python实现基于米勒拉宾素性检测算法最快的超大数素性检测

python, 快速幂模, 最快的超大数素性检测, 米勒拉宾素性检测算法, 蒙哥马利算法, 超大数素性检测

0x00 米勒拉宾素性检测算法伪代码

米勒拉宾素性检测算法的伪代码：

Input #1: n > 3, an odd integer to be tested for primality;
Input #2: k, a parameter that determines the accuracy of the test
Output: composite if n is composite, otherwise probably prime

write n − 1 = (2 ^ r) * d with d odd by factoring powers of 2 from n − 1
WitnessLoop: repeat k times:
   pick a random integer a in the range [2, n − 2]
   x ← a^d mod n
   if x = 1 or x = n − 1 then
      continue WitnessLoop
   repeat r − 1 times:
      x ← x^2 mod n
      if x = 1 then
         return composite
      if x = n − 1 then
         continue WitnessLoop
   return composite
return probably prime

0x01 实现

一开始针对这个米勒拉宾素性检测算法里的幂模运算，我以为用蒙哥马利算法－快速幂模要比pow()函数更快，结果还是内置的pow()更快：

1
2
3
4
5
6
7
8
9
# x = a^d n n
def montgomery(a, d, n):
    x = 1
    while d:
        if (d & 1):
            x = (x * a) % n
        d >>= 1  
        a = (a * a) % n
    return x

上面的伪代码中a的选取是可以优化的，可以看到下面代码有实现，原理看一下wikipedia：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
def select_a_list(n):
    if n < 2047:
        return [2,]
    elif n < 1373653:
        return [2, 3]
    elif n < 9080191:
        return [31, 73]
    elif n < 25326001:
        return [2, 3, 5]
    elif n < 3215031751:
        return [2, 3, 7]
    elif n < 4759123141:
        return [2, 7, 61]
    elif n < 1122004669633:
        return [2, 13, 23, 1662803]
    elif n < 2152302898747:
        return [2, 3, 5, 7, 11]
    elif n < 3474749660383:
        return [2, 3, 5, 7, 11, 13]
    elif n < 341550071728321:
        return [2, 3, 5, 7, 11, 13, 17]
    elif n < 3825123056546413051:
        return [2, 3, 5, 7, 11, 13, 17, 19, 23]
    elif n < 318665857834031151167461:
        return [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
    elif n < 3317044064679887385961981:
        return [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]
    else:
        return [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]

0x02 测试

用这个方法检测素数，100位10进制数平均0.0003，1000位10进制数平均0.09秒，10000位10进制数平均60秒。

0x03 完整代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
#!/usr/bin/env python
# coding=utf-8
# author:admin[@hackfun.org]
# license:GPL v3
# blog:hackfun.org
def select_a_list(n):
    if n < 2047:
        return [2,]
    elif n < 1373653:
        return [2, 3]
    elif n < 9080191:
        return [31, 73]
    elif n < 25326001:
        return [2, 3, 5]
    elif n < 3215031751:
        return [2, 3, 7]
    elif n < 4759123141:
        return [2, 7, 61]
    elif n < 1122004669633:
        return [2, 13, 23, 1662803]
    elif n < 2152302898747:
        return [2, 3, 5, 7, 11]
    elif n < 3474749660383:
        return [2, 3, 5, 7, 11, 13]
    elif n < 341550071728321:
        return [2, 3, 5, 7, 11, 13, 17]
    elif n < 3825123056546413051:
        return [2, 3, 5, 7, 11, 13, 17, 19, 23]
    elif n < 318665857834031151167461:
        return [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
    elif n < 3317044064679887385961981:
        return [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]
    else:
        return [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]
def calc_r_d(n):
    n -= 1
    r = 0
    d = n
    while 1:
        if n % 2 != 0:
            break
        r += 1
        d = n / 2
        n /= 2
    return r, d
def miller_rabin_primality_test(n):
    """use miller rabin primality test to judge n whether prime"""
    if n < 0:
        n = -n
    if n < 2:
        return False
    if n == 2 or n == 3:
        return True
    r, d = calc_r_d(n)
    a_list = select_a_list(n)
    for k in xrange(len(a_list)):
        a = random.sample(a_list, 1)[0] # Select one non-repeating random number from the list at a time
        a_list.remove(a)
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        for t in xrange(r - 1):
            x = pow(x, 2, n)
            if x ==1:
                return False
            if x == n - 1:
                break
        else:
            return False
    return True