从前有一种数,他们很孤独,只有1这唯一的朋友,他们彼此间也找不到交集,他们永远都在那里躺着。对,他们就是素数,我们应该背过一个叫做素数表的东西,唔。。只是100以内的。
让我们回顾一下: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。所以如何找到这些坏家伙呢,这个问题一直困扰着我们,欧拉,费马,高斯等等超级大神也都对素数抱有极大的兴趣。
下面让我们来看看几种素数的判定方法:
最简单的素数判定
既然素数是只能被1和自身整除的正整数,那么最直接的手段当然是一个个去试咯,如果能被其他数整除,就不是素数了。测试的时候只需要测试到
floor(n√)
即可,因为大于
floor(n√)
的要出现早就已经出现了
于是我们自然会想到这样的代码:
def Isprime(n):
return 0 not in [n%i for i in range(2,int(n**0.5)+1)]
可以优化一点点:不能整除2,那肯定2 的倍数也就不用考虑了。所以只考虑奇数和2就行。(偶数和1的与为0)
def Isprime(n):
if n==2: return True
elif n<2 or n&1==0: return False
return 0 not in (n%i for i in range(3,int(n**0.5)+1,2))
然而这样只是减小了一半,并没有实质性的效果,对于需要测试的n来说,复杂度依然是
O(n√)
级别。有没有一种常数级别的测试方法呢?
答案是显然的,超级大神们给出了一系列的定理。让我们来看看:
费马小定理
费马小定理
a是一个整数,p是一个质数,那么 ap−a 是 p 的倍数,可以表示为
如果a不是p的倍数,这个定理也可以写成ap≡a(modp)
ap−1≡1(modp)证明也很简单:
考虑 [1:p−1]∗a 这 p−1 个数,这个向量除以 p 的余数组成的集合一定是[1:p−1] ,这是因为对于任两个相异 k∗a 而言 (k=1,2,3....(p−1)) ,其差不是p的倍数,且任一个 k∗a 亦不为p的倍数(所以余数不为0)。因此:
1⋅2⋅3⋅⋯⋅(p−1)≡(1⋅a)⋅(2⋅a)⋅⋯⋅((p−1)⋅a)(modp)两边约去 (p−1)! 即得:
ap−1≡1(modp)
有着如此强大的定理撑腰,我们就可以考虑新的算法了。
Miller-Rabin算法
新算法也很简单,随便找几个数
a
,看
这种算法在n比较大的时候有极大的优势,因为运算次数并不会因为n变大而变大多少,相比之前的
O(n√)
级别的增长,这个可以看作是常数时间。然而事实果真如此吗?
我们考虑一下算法中提到的IsP测试,它需要进行一次指数运算和一次求模。当n很大的时候,
ap−1
会溢出,随便举个例子吧,今天是12月12日,如果我们想测试1213的素性(正好是个素数),则我们很有可能取到
12121213
这个数,唔。。我们可以大概估算一下这个数有多大,用计算机存这个数的话应该需要
log121212132=1213log12122=12424.97
位才能存下。唔。。是不是太大了。。
于是问题又化成了如何快速求
ab%m
,我们需要这样一个定理:
积模分解定理:
(a∗b)%m=((a%m)∗(b%m))%m
写成代码就是:
def ExpMod(a,b,m):#return a**b%m
k=1
a%=m
while b!=1:
if b&1!=0:
k=(k*a)%m
a=(a*a)%m
b>>=1
return (a*k)%m
解决了这个问题,我们的Miller-Rabin素数判定算法就可以直接写了:
测试次数取50,当然越大越好,但是一般来说已经足够。
from random import randint
def Miller_Rabin(n,s=50):
'''s means test s times,the higher the better but slower as well'''
if n==2: return True
elif n<2 or n&1==0: return False
for i in range(s):
a=randint(1,n-1)
if ExpMod(a,n-1,n)!=1:
return False
return True
复杂度分析:由于Expmod算法是
O(log2n)
的,而测试一共进行了s次,为常数时间,故总的复杂度为
O(log2n)
,鉴于有常数因子的存在,故在数量比较小的时候普通的Isprime算法要比较好。
说了这么多,上个测试吧:
from time import time
def test(n,t):
t1=time()
for i in range(t):
lis1=Isprime(n)
t2=time()
for i in range(t):
lis2=Miller_Rabin(n)
t3=time()
print('n=',n,'\ntimes=',t,'\nIsprime:',lis1,'t=',t2-t1,'\nMiller:',lis2,'t=',t3-t2)
test(3215031751,100)
test(3215031751111,10)
结果如下:
n= 3215031751
times= 100
Isprime: False t= 0.002001047134399414
Miller: True t= 0.09900498390197754n= 3215031751111
times= 10
Isprime: True t= 1.7500998973846436
Miller: True t= 0.014001131057739258
咦~发现了什么没有:3215031751 竟然用两种方法测试出来结果不一样。这是为什么呢?
我们知道,Isprime算法一定是准确无误的,而Miller-Rabin是一个概率算法,而恰恰就有这样的数C,C本身是合数,但是依然能通过测试,这种数被称为Carmichael数。然而,在Miller算法中已经通过随机选取a来避免误判了,然而还是会有这种漏网之鱼,这种时候我们选取较大的s,即加大测试次数即可。
另外我们发现这个漏网之鱼相对来说是比较小的,在这种情况下,Isprime算法要快一些,因为他只要测试到合数的第一个质因子就结束了,然而对于合数Miller-Rabin算法会一直测试所有的s次, 所以在n比较小的时候,应该选取Isprime算法。
但是当n比较大,而且是质数的时候,比如下面的3215031751111,这是一个质数,Isprime则需要计算到
floor(n∗∗0.5)+1=1793051
次才能结束,Miller算法依然只有50多次,相比之下,时间也是相差了上百倍。则对于比较大的n,Miller-Rabin算法是可信的。
已经证明:Miller-Rabin算法一定能确定素数,而其误把合数判为素数的概率为 2−s 。
生成素数表
前面的讨论中给出了两种常见的素数判定方法,那么如果我们想得到n素数表,该如何操作呢。有人会说,从2到n每个奇数挨着判断呗,这是一个办法,然而,考虑到对每个数都会进行素数判断,实在是略麻烦。
不如反过来考虑,我们可以在数轴上挖坑,先把2的倍数全挖掉,这时,我们剩下的最小的数是3,那么3就是一个质数;继续刚才的操作,不过这时我们要开始挖3的倍数,这时我们剩下了5,唔。。看来5也是一个质数,那继续挖5呢,剩下了7。。。所以,这么挖下去,每次取最小的剩下的那个不就是下一个素数了么。就这么一直挖下去,每次挖都用新得到的素数,我们就可以得到小于n的所有的素数了。
废话少说,看图:
这个算法叫Sieve Eratosthenes.就上C#代码吧。。返回一个迭代器。以前写的懒得改了。。=。=有时间再写个Python版本。。
static IEnumerable<long> EraSieve(long n)
{
bool[] isPrime=new bool[n+1];
for(int i=3;i<=n;i+=2)
isPrime[i]=true;
isPrime[2]=true;
for(int i=3;i<=Math.Sqrt(n);i+=2)
{
if(isPrime[i]==true)
{
for(int j=i*i;j<=n;j+=2*i)
isPrime[j]=false;
}
}
for(int i=1;i<=n;i++){
if(isPrime[i])
yield return i;
}
}