# 数论函数

## 最大公约数

### 程序设计



In [3]:

def gcd(a, b):
# returns GCD with Euclid's Algorithm
while a != 0:
a, b = b % a, a
return b



## 模逆

### 程序设计



In [4]:

def modInverse(a, m):
# a and m must be relatively prime
if gcd(a, m) != 1:
return None
# Extended Euclidean Algorithm
u1, u2, u3 = 1, 0, a
v1, v2, v3 = 0, 1, m
while v3 != 0:
q = u3 // v3 # // is the integer division operator
v1, v2, v3, u1, u2, u3 = \
(u1 - q * v1), (u2 - q * v2), (u3 - q * v3), v1, v2, v3
return u1 % m



## Miller Rabin 素性检测

### 程序设计



In [5]:

import random

def millerRabin(n):
# Returns True if prime
# n - 1 = 2^k * m, m is odd
m = n - 1
k = 0
while m % 2 == 0:
# floor division
m //= 2
k += 1
# random number a (1 <= a <= n-1)
a = random.randrange(1, n)
# b = a^m mod n
# pow() with 3 arguments is faster than 2
b = pow(a, m, n)

if b % n == 1:
return True

for i in range(0, k):
if (b + 1) % n == 0:
return True
else:
b = pow(b, 2, n)

return False



## 复合素数检测函数

### 算法

1. 通过素数筛法(classicPrime.py)生成 1000 以内的所有素数列表。
在内部以列表的形式存储

2. 若待检测数在素数表内，就是素数

3. 若素数表内素数为待检测数因子，即为合数

4. 上述两点判断皆失效时，调用 ***Miller Rabin 素性检测***



### 程序设计

• Eratosthenes 素数筛函数


In [6]:

%%file classicPrime.py

import math
# Determine Prime with the help of
# the Sieve of Eratosthenes algorithm
def primeSieve(sieveSize):

# Returns a list of prime numbers
sieve = [True] * sieveSize
sieve[0] = sieve[1] = False

# create the sieve
for i in range(2, int(math.sqrt(sieveSize)) + 1):
multiplier = 2
while i * multiplier < sieveSize:
sieve[i * multiplier] = False
multiplier += 1

# create list
primelist = []
for i in range(sieveSize):
if sieve[i] == True:
primelist.append(i)

return primelist




Overwriting classicPrime.py


• 复合素数检测函数


In [26]:

def isPrime(num):

if (num < 2):
return False

lowPrimes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, \

41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103,\

107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, \

173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, \

239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, \

311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, \

383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, \

457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, \

541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, \

613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, \

683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, \

769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, \

857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, \

941, 947, 953, 967, 971, 977, 983, 991, 997]

if num in lowPrimes:
return True

for prime in lowPrimes:
if (num % prime == 0):
return False

return millerRabin(num)


• 概率型大素数生成函数

接受传入的素数位数（二进制），默认1024。

生成随机数，调用素数检测，检测成功则返回。



In [8]:

def generateLargePrime(keySize=1024):
while True:
num = random.randrange(2**(keySize-1), 2**(keySize))
if isPrime(num):
return num



# RSA 参数生成

• 按照秘钥长度生成参数

按照教材密码体制 5.1 实现

默认使用私钥加密， 公钥解密。因此，在生成密钥时，保存私钥的一方同时保存 $p , q$ 两个奇素数，为使用中国剩余定理进行加速实现提供方便。



In [9]:

import random, time

def parameterGenerator(keySize=1024):
time0 = timenow = time.time()
# Step 1 : Generate two random large primes
p = generateLargePrime(keySize)
q = generateLargePrime(keySize)
n = p * q
# phi(n) -> phin, Euclidean Phi Function
phin = (p - 1) * (q - 1)

# Step 2 :
# Generate a number relatively prime to phi(n)
timenow = time.time()
while True:
# Keep trying random numbers
e = random.randrange(2 ** (keySize -1 ), 2 ** keySize)
# if e and phi(n) are relatively prime,
# stop trying new e, break the loop
if gcd(e, phin) == 1:
break

# Step 3 : Calculate d, the mod inverse of e
timenow = time.time()
d = modInverse(e, phin)
print("""(n, p, q, e, d)""")
return (n, p, q, e, d)



# RSA 原始实现

加解密内容均为整数


• 加密函数


In [10]:

def rsaPrimitiveEncryption(message, n, e):
encryptedMessage = pow(message, e, n)
return encryptedMessage


• 解密函数


In [11]:

def rsaPrimitiveDecryption(message, n, d):
decryptedMessage = pow(message, d, n)
return decryptedMessage



# RSA 快速实现

## 使用中国剩余定理加速解密

### 中国剩余定理数学函数

• 问题描述

• 需要求解如下方程：

$$x \equiv a \pmod m$$ $$x \equiv b \pmod n$$ $$(m, n) = 1$$

• 需要给出如下结果：

$$x \equiv K \pmod M$$ $$M = a * b$$

• IO

- 使用字典数据结构作为传入参数:

modSys = {m: a, n: b}

- 使用二维数组作为传出参数：

result = (M, K)


• 异常输入处理（未处理）

• 除数为零

• 除数不互素

• 测试用例

• 需要求解如下方程：

$$x \equiv 1 \pmod 5$$ $$x \equiv 3 \pmod 7$$

• 需要给出如下结果：

$$x \equiv 31 \pmod{35}$$

• 算法

• 程序设计


In [12]:

def crt(modSys = {5: 1, 7: 3}):
# build divisors list and divisors producat
M, divisors = 1, []
for divisor in modSys:
M *= divisor
divisors.append(divisor)
# start solving by adding divisors
K = 0;
for i in divisors:
K += modSys[i]*int(M/i)*modInverse(a=int (M/i), m=i)
K = int(K % M)
result = (M, K)
return result



### 中国剩余定理实现的 RSA 加速解密函数

• 原始实现

求解同余式：

$$P = C^{d} \pmod{n}$$

• 中国剩余定理

求解同余式组：

$$P = C^{d} \pmod{p}$$

$$P = C^{d} \pmod{q}$$



In [13]:

def rsaCrtDecryption(message, n, p, q, d):
r1 = pow(message, d, p)
r2 = pow(message, d, q)
(M, K) = crt({p: r1, q: r2})
if M == n:
decryptedMessage = K
return decryptedMessage
else:
print("Fatal Error")
return None



## 使用模重复平方加速解密

### 模重复平方数学函数

• 问题描述

当$k$ 很大时，需要计算：

$$a^k \pmod{m}$$

• 算法实现

按照 A Friendly Introduction to Number Theory 算法16.1 实现

1. 将 $k$ 表示为 $2$ 的幂次和（ $k$ 的二进制展开）

2. 使用模重复平方法制作模 $m$ 的 $a$ 的幂次表

3. 对上述幂次进行累乘

• IO

• 输入

a, k, m

• 输出

remainder

• 测试用例

• 求解

$$7^{327} \pmod{853}$$

• 预计结果

$$7^{327} \equiv 286 \pmod{853}$$



In [14]:

from functools import reduce
def successiveSquaring(a, k, m):
# binaryExpansion
binExpList = []
for(offset, bit) in enumerate(bin(k).replace('0b','')):
binExpList.append(int(bit, 2))
binExpList.reverse()

# successiveSquaring
powersOfAModuloMList = []
for bit in range(len(binExpList)):
if bit == 0:
powersOfAModuloMList\
.append(pow(a, pow(2, bit), m))
else:
powersOfAModuloMList\
.append(pow(powersOfAModuloMList[bit - 1], 2, m))
# product
l = \
list(map(lambda bit,powersOfAModuloM: \
bit*powersOfAModuloM,\
binExpList, powersOfAModuloMList))
while 0 in l:
l.remove(0)
return reduce(lambda x,y:x*y%m, l)



#### 模重复平方实现的 RSA 加速加解密函数

• 加密函数


In [15]:

def rsaSecSqrEncryption(message, n, e):
encryptedMessage = successiveSquaring(message, e, n)
return encryptedMessage


• 解密函数


In [16]:

def rsaSecSqrDecryption(message, n, d):
decryptedMessage = successiveSquaring(message, d, n)
return decryptedMessage



## 使用蒙哥马利算法加速解密

### 蒙哥马利算法数学函数

• 问题描述

当$e$ 很大时，需要计算：

$$a^e \pmod{m}$$

• IO

• 输入

a, e, m

• 输出

remainder

• 算法实现

按照 Cryptography in C and C++ $\S 6.4$ 实现

• 测试用例

• 求解

$$7^{327} \pmod{853}$$

• 预计结果

$$7^{327} \equiv 286 \pmod{853}$$

• 算法说明

• $n' \equiv - n^{-1} \pmod{r}$

• $r' \equiv r^{-1} \pmod{n}$

• 程序设计


In [17]:

%%file montgomery.cpp

#include "gmp.h"
#define TESTCOUNT 500
#define TESTBITS 1024
#define MONTBITS (8*sizeof(mp_limb_t))
#define MONT_MAX 64

gmp_randstate_t state;

void MontMulti(mpz_t T,
const mpz_t x,
const mpz_t y,
const  mpz_t N,
const mp_limb_t N_1)
{
unsigned int i;
mp_limb_t num, carry, res[MONT_MAX] = { 0 };
mp_limb_t *temp,t[MONT_MAX] = { 0 };

if (x->_mp_size > y->_mp_size)
mpn_mul(t, x->_mp_d, x->_mp_size,y->_mp_d,y->_mp_size);
else
mpn_mul(t, y->_mp_d, y->_mp_size, x->_mp_d, x->_mp_size);

temp = t;
for (i = 0; i < N->_mp_size; i++){
num = temp[0]*N_1;
res[i]=mpn_addmul_1(temp, N->_mp_d,N->_mp_size,num);
temp++;
}

carry = mpn_add_n(temp, temp, res, N->_mp_size);
if (carry != 0 || mpn_cmp(temp, N->_mp_d, N->_mp_size) >= 0)
mpn_sub_n(temp, temp, N->_mp_d, N->_mp_size);

mpz_import(T, N->_mp_size, -1, sizeof(mp_limb_t), 0, 0,temp);
}
void Mont_Exp(mpz_t rop,
const mpz_t base,
const mpz_t exp,
const mpz_t N)
{
mp_limb_t N_1;
mpz_t K, P, R, temp, N_inv, b;
mp_bitcnt_t index;
unsigned int bitnum = mpz_sizeinbase(exp, 2);
unsigned int rbit = N->_mp_size*MONTBITS;
mpz_inits(K, P, R, temp, N_inv, b, NULL);

mpz_setbit(b, MONTBITS);
mpz_invert(N_inv, N, b);
mpz_sub(N_inv, b, N_inv);
N_1 = *(N_inv->_mp_d);
mpz_set_ui(temp, 1);
mpz_setbit(K, 2 * rbit);
mpz_mod(K, K, N);

MontMulti(P, K, base, N, N_1);
MontMulti(R, K, temp, N, N_1);
for (index = 0; index < bitnum; index++) {
if (mpz_tstbit(exp, index) == 1)
MontMulti(R, R, P, N, N_1);
MontMulti(P, P, P, N, N_1);
}
MontMulti(rop, temp, R, N, N_1);

mpz_clears(K, P, R, temp, N_inv, b, NULL);
}

int main(int argc, char* argv[])
{
gmp_randinit_lc_2exp_size(state, 128);
gmp_randseed_ui(state, (unsigned long)time(NULL));
int i;
time_t t1, t2;
mpz_t base, exp, mod, r1, r2;
mpz_inits(base, exp, mod, r1, r2, NULL);

mpz_rrandomb(exp, state, TESTBITS);
mpz_rrandomb(mod, state, TESTBITS);
mpz_rrandomb(base, state, TESTBITS);
mpz_setbit(mod, 0);
mpz_mod(base, base, mod);

for (i = 0; i < TESTCOUNT; i++) {
Mont_Exp(r1, base, exp, mod);
}
mpz_clears(base, exp, mod, r1, r2, NULL);
return 0;
}




Overwriting montgomery.cpp



# RSA 原始及快速实现时间比对测试

1 原始实现 Python 15.2s 500
2 中国剩余定理 Python 9.75s 500
3 模重复平方 Python 5.51s 500
4 蒙哥马利算法 C++ 0.738s 500

• 使用三种加速算法显著加速了RSA算法的解密速度。

• 而使用 C++ 及 GMP 库实现的算法弥补了Python 脚本化运行的效率短板，将速度提升了一个数量级。

## 原始实现



In [18]:

(n, p, q, e, d) = parameterGenerator()
plain = random.randint(2**10, 2**11)
cipher = rsaPrimitiveEncryption(e=e,message=plain, n=n)




(n, p, q, e, d)




In [19]:

%%time
for i in range(500):
rsaPrimitiveDecryption(d=d,message=cipher, n=n)




CPU times: user 15.2 s, sys: 19.9 ms, total: 15.2 s
Wall time: 15.2 s



## 中国剩余定理



In [20]:

(n, p, q, e, d) = parameterGenerator()
plain = random.randint(2**10, 2**11)
cipher = rsaPrimitiveEncryption(e=e,message=plain, n=n)




(n, p, q, e, d)




In [21]:

%%time
for i in range(500):
rsaCrtDecryption(d=d, message=cipher, n=n,p=p,q=q)




CPU times: user 9.73 s, sys: 14.2 ms, total: 9.74 s
Wall time: 9.75 s



## 模重复平方



In [22]:

(n, p, q, e, d) = parameterGenerator()
plain = random.randint(2**10, 2**11)
cipher = rsaPrimitiveEncryption(e=e,message=plain, n=n)




(n, p, q, e, d)




In [23]:

%%time
for i in range(500):
rsaSecSqrDecryption(d=d, message=cipher, n=n)




CPU times: user 5.5 s, sys: 8.52 ms, total: 5.51 s
Wall time: 5.51 s



## 蒙哥马利算法



In [24]:

!g++ ./montgomery.cpp -l gmp




In [25]:

%%time
!./a.out




CPU times: user 9.68 ms, sys: 7.88 ms, total: 17.6 ms
Wall time: 738 ms