数论函数

最大公约数

算法

使用欧几里得算法

程序设计



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