About素数

　　从前有一种数，他们很孤独，只有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。所以如何找到这些坏家伙呢，这个问题一直困扰着我们，欧拉，费马，高斯等等超级大神也都对素数抱有极大的兴趣。
　　下面让我们来看看几种素数的判定方法：

最简单的素数判定

　　既然素数是只能被１和自身整除的正整数，那么最直接的手段当然是一个个去试咯，如果能被其他数整除，就不是素数了。测试的时候只需要测试到 $floor(\sqrt{n})$ 即可，因为大于 $floor(\sqrt{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(\sqrt{n})$级别。有没有一种常数级别的测试方法呢？
　　答案是显然的，超级大神们给出了一系列的定理。让我们来看看：

费马小定理

费马小定理

a是一个整数，p是一个质数，那么 $a^{p}-a$ 是 $p$ 的倍数，可以表示为

$$a^{p}\equiv a{\pmod {p}}$$ 如果a不是p的倍数，这个定理也可以写成

$$a^{p-1}\equiv 1{\pmod {p}}$$ 证明也很简单：
　　考虑 $[1:p-1]*a$ 这 $p-1$ 个数，这个向量除以 $p$ 的余数组成的集合一定是 $[1:p-1]$ ，这是因为对于任两个相异 $k*a$ 而言 $(k=1,2,3....(p-1))$，其差不是p的倍数，且任一个 $k*a$ 亦不为p的倍数（所以余数不为0）。因此：
$$1\cdot 2\cdot 3\cdot \dots \cdot (p-1)\equiv (1\cdot a)\cdot (2\cdot a)\cdot \dots \cdot ((p-1)\cdot a){\pmod {p}}$$ 两边约去 $(p-1)!$ 即得：
$$a^{p-1}\equiv 1{\pmod {p}}$$

　　有着如此强大的定理撑腰，我们就可以考虑新的算法了。

Miller-Rabin算法

　　新算法也很简单，随便找几个数$a$，看 $a^{p-1}\% p$ 是否等于1呗。我们记这个判断为IsP。于是问题就转换成了找一些数 $a\in[1:p-1]$，如果这些数均通过了IsP测试，则我们有极大的把握认为p为质数。这样一来，我们只需要计算固定次数的指数和模运算，便可以有极大的把握判定一个数的素性。这个算法称为Miller-Rabin算法。
　　这种算法在n比较大的时候有极大的优势，因为运算次数并不会因为n变大而变大多少，相比之前的$O(\sqrt{n})$级别的增长，这个可以看作是常数时间。然而事实果真如此吗？
　　我们考虑一下算法中提到的IsP测试，它需要进行一次指数运算和一次求模。当n很大的时候，$a^{p-1}$ 会溢出，随便举个例子吧，今天是12月12日，如果我们想测试1213的素性（正好是个素数），则我们很有可能取到 $1212^{1213}$ 这个数，唔。。我们可以大概估算一下这个数有多大，用计算机存这个数的话应该需要 $\log_{2}^{1212^{1213}}=1213\log_{2}^{1212}=12424.97$ 位才能存下。唔。。是不是太大了。。
　　于是问题又化成了如何快速求 $a^{b}\%{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(\log_{2}{n})$ 的，而测试一共进行了s次，为常数时间，故总的复杂度为$O(\log_{2}{n})$，鉴于有常数因子的存在，故在数量比较小的时候普通的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.09900498390197754

n= 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}$ 。

生成素数表

　　前面的讨论中给出了两种常见的素数判定方法，那么如果我们想得到ｎ素数表，该如何操作呢。有人会说，从２到ｎ每个奇数挨着判断呗，这是一个办法，然而，考虑到对每个数都会进行素数判断，实在是略麻烦。
　　不如反过来考虑，我们可以在数轴上挖坑，先把２的倍数全挖掉，这时，我们剩下的最小的数是３，那么３就是一个质数；继续刚才的操作，不过这时我们要开始挖３的倍数，这时我们剩下了５，唔。。看来５也是一个质数，那继续挖５呢，剩下了７。。。所以，这么挖下去，每次取最小的剩下的那个不就是下一个素数了么。就这么一直挖下去，每次挖都用新得到的素数，我们就可以得到小于ｎ的所有的素数了。
　　废话少说，看图：
　　
　　这个算法叫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;
            }
        }