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