素数的Miller_Rabbin测试

【引理】如果p是一个素数的话，那么对任意一个小于p的正整数a，a, 2a, 3a, ..., (p-1)a除以p的余数正好是一个1到p-1的排列。（例如，5是素数，3, 6, 9, 12除以5的余数分别为3, 1, 4, 2，正好就是1到4这四个数。）

  【证明】（ 反证法）假如结论不成立的话，那么就是说有两个小于p的正整数m和n使得na和ma除以p的余数相同。不妨假设n>m，则p可以整除a(n-m)。但p是素数，那么a和n-m中至少有一个含有因子p。这显然是不可能的，因为a和n-m都比p小。

【Fermat小定理（Fermat's Little Theorem）】如果p是素数，a是小于p的正整数，那么a^(p-1) mod p ≡ 1。

  【证明】用同余式表述，可以证明：(p-1)! ≡ a * 2a * 3a * ... * (p-1)a (mod p)，也即：(p-1)! ≡ (p-1)! * a^(p-1) (mod p)。

两边同时除以(p-1)!，就得到了我们的最终结论：1 ≡ a^(p-1) (mod p)

 

---

 

谈到Fermat小定理，数学历史上有很多误解。很长一段时间里，人们都认为Fermat小定理的逆命题是正确的，并且有人亲自验证了a=2, p<300的所有情况。1819年有人发现了Fermat小定理逆命题的第一个反例：虽然2的340次方除以341余1，但341=11*31。后来，人们又发现了561, 645, 1105等数都表明a=2时Fermat小定理的逆命题不成立。虽然这样的数不多，但不能忽视它们的存在。于是，人们把所有能整除2^(n-1)-1的合数n叫做伪素数(Pseudoprime数)。

不满足2^(n-1) mod n = 1的n一定不是素数；如果满足的话则多半是素数。这样，一个比试除法效率更高的素性判断方法出现了：制作一张伪素数表，记录某个范围内的所有伪素数，那么所有满足2^(n-1) mod n = 1且不在伪素数表中的n就是素数。之所以这种方法更快，是因为我们可以使用二分法快速计算2^(n-1)  mod n 的值，这在计算机的帮助下变得非常容易；在计算机中也可以用二分查找有序数列、Hash表开散列、构建Trie树等方法使得查找伪素数表效率更高。

由于伪素数的存在使得我们需要改进素数检测算法。最简单的想法就是，我们刚才只考虑了a=2的情况。对于式子a^(n-1) mod n，取不同的a可能导致不同的结果。一个合数可能在a=2时通过了测试，但a=3时的计算结果却排除了素数的可能。于是，人们扩展了伪素数的定义，称满足a^(n-1) mod n = 1的合数n叫做以a为底的伪素数(pseudoprime to base a)。前10亿个自然数中同时以2和3为底的伪素数只有1272个，这个数目不到刚才的1/4。这告诉我们如果同时验证a=2和a=3两种情况，算法出错的概率降到了0.000025。容易想到，选择用来测试的a越多，算法越准确。通常我们的做法是，随机选择若干个小于待测数的正整数作为底数a进行若干次测试，只要有一次没有通过测试就立即把这个数扔回合数的世界。这就是Fermat素性测试。

然而并不能仅通过选取所有小于n的基数a来消除素数测试中的出错机会，因为总存在这样的合数，它可以通过所有的测试。Carmichael第一个发现这样极端的伪素数，因此这类数被称作Carmichael数。Carmichael数的存在说明，我们还需要继续加强素性判断的算法。

Miller和Rabin两个人的工作让Fermat素性测试迈出了革命性的一步，建立了Miller-Rabin素性测试算法。新的测试基于下面的定理：如果p是素数，x是小于p的正整数，且x² mod p = 1，那么要么x=1，要么x=p-1。这是显然的，因为x² mod p = 1相当于p能整除x²-1，也即p能整除(x+1)(x-1)。由于p是素数，那么只可能是x-1能被p整除(此时x=1)或x+1能被p整除(此时x=p-1)。

这就是Miller-Rabin素性测试的方法：不断地提取指数n-1中的因子2，把n-1表示成d*2^r（其中d是一个奇数）。那么我们需要计算的东西就变成了a的d*2^r次方除以n的余数。于是，a^(d * 2^r-1)要么等于1，要么等于n-1。如果a^(d *2^r-1)等于1，定理继续适用于a^(d * 2^r-2)，这样不断开方开下去，直到对于某个i满足a^(d * 2ⁱ) mod n = n-1或者最后指数中的2用完了得到的a^dmod n=1或n-1。

这样，Fermat小定理加强为如下形式：尽可能提取因子2，把n-1表示成d*2^r，如果n是一个素数，那么或者a^d mod n=1，或者存在某个i使得a^(d*2ⁱ) mod n=n-1 ( 0<=i<r ) （注意i可以等于0，这就把a^d mod n=n-1的情况统一到后面去了）.

Miller-Rabin素性测试同样是不确定算法，我们把可以通过以a为底的Miller-Rabin测试的合数称作以a为底的强伪素数(strong pseudoprime)。第一个以2为底的强伪素数为2047。第一个以2和3为底的强伪素数则大到1 373 653。

Miller-Rabin算法的简单实现（C语言）如下：

long exp_mod(long a,long b,long n){
    long r=1;
    while(b){
        if(b&1)r=(r*a)%n;
        a=(a*a)%n;
        b>>=1;        
    }
    return r;
}
bool Miller_Rabbin(long n){
    long a,s=(n-1)>>2,i=1,t=n-1;
    do{
        a=rand()%(t)+1;
        if(exp_mod(a,t,n)!=1)return false;
    }while(i++<s);
    return true;
}