素性测试

2014.12.16 17:46 Tue| 5 visits oi_2015| 2015_算法笔记| Text

素性测试出错的概率相当于你被闪电击中的概率。 ——敬部长一杯酒

那么如果我要是跑一晚上素性测试而且我没有被闪电击中,那么是不是说素性探测没有出错? ——某沙茶

米勒-拉宾素性检验

米勒-拉宾素性检验是一种素数判定法则,利用随机化算法判断一个数是合数还是可能是素数。卡内基梅隆大学的计算机系教授 Gary Lee Miller 首先提出了基于广义黎曼猜想的确定性算法,由于广义黎曼猜想并没有被证明,其后由以色列耶路撒冷希伯来大学的Michael O. Rabin教授作出修改,提出了不依赖于该假设的随机化算法。

要测试 $N$ 是否为素数,首先将 $N-1$ 分解为 $2^t u$ 。在每次测试开始时,先随机选一个介于 $[1, N-1]$ 的整数 $a$ ,之后如果对所有的 $r \in [0, t-1]$ ,若 $a^u \mod N \neq 1$ 且 $a^{2^{r}u} \mod N \neq -1$ ,则 $N$ 是合数。否则,$N$ 有 $3/4$ 的概率为素数。

摘自中文wikipedia

如果MILLER-RABIN返回PRIME,则它仍有一种很小的可能性会产生错误。然而,出错的可能性并不依赖于n,对该过程也不存在坏的输入。相反,它取决于次数和在选取基值a时“抽签的运气”。因此从总的原则上,对随机选取的整数n,其出错率应该是很小的。

摘自算法导论

算法实测

下面的代码计算出了形如2333***333的素数(3的个数在1000以内)。结果如下(数据中的数字表示3的个数):

1 may be a prime
2 may be a prime
3 may be a prime
4 may be a prime
10 may be a prime
16 may be a prime
22 may be a prime
53 may be a prime
91 may be a prime
94 may be a prime
106 may be a prime
138 may be a prime
210 may be a prime
282 may be a prime
522 may be a prime
597 may be a prime

测试在一次测试中,100以内的素数的识别正确性,证明一次测试的错误大约为3个,多为9、15、65等素数出现错误。在两次测试中仅出现一次在15处出现错误。可见素数测试的正确率的确较高。

代码实现

import math
import random

def power(a,b,p):
    re = 1
    while(b > 0):
        if(b & 1):
            re = re*a%p
        a = a*a%p
        b = b>>1
    return re

def witness(a,n):
    u=n-1
    t=0
    while(u&1 == 0):
        u = u>>1
        t = t+1
    x={}
    x[0] = power(a,u,n)
    for i in range(1,t+1):
        x[i] = x[i-1]**2 % n
        if(x[i]==1 and x[i-1]!=1 and x[i-1]!=n-1):
            return 1
    if(x[i]!=1):
        return 1
    return 0

def isprime(n,s):
    if(n==2): return 1
    if(n&1 == 0): return 0
    for i in range(1,s+1):
        if(witness(random.randint(1,n-1),n) == 1):
            return 0
    return 1

n=2
for i in range(1,1001):
    n = n*10+3
    if(isprime(n,100) == 1):
        print(n,"may be a prime")