标签：style blog http color os strong io 数据
首先附上matrix67大神的讲解：
-----------------------------------------------------------------------------------------------------------------------------------------------------
Miller和Rabin两个人的工作让Fermat素性测试迈出了革命性的一步，建立了传说中的Miller-Rabin素性测试算法。新的测试基于下面的定理：如果p是素数，x是小于p的正整数，且x^2 mod p = 1，那么要么x=1，要么x=p-1。这是显然的，因为x^2 mod p = 1相当于p能整除x^2-1，也即p能整除(x+1)(x-1)。由于p是素数，那么只可能是x-1能被p整除(此时x=1)或x+1能被p整除(此时x=p-1)。
我们下面来演示一下上面的定理如何应用在Fermat素性测试上。前面说过341可以通过以2为底的Fermat测试，因为2^340 mod 341=1。如果341真是素数的话，那么2^170 mod 341只可能是1或340；当算得2^170 mod 341确实等于1时，我们可以继续查看2^85除以341的结果。我们发现，2^85 mod 341=32，这一结果摘掉了341头上的素数皇冠，面具后面真实的嘴脸显现了出来，想假扮素数和我的素MM交往的企图暴露了出来。
这就是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^i) mod n = n-1或者最后指数中的2用完了得到的a^d mod n=1或n-1。这样，Fermat小定理加强为如下形式：
尽可能提取因子2，把n-1表示成d*2^r，如果n是一个素数，那么或者a^d mod n=1，或者存在某个i使得a^(d*2^i) 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算法的代码也非常简单：计算d和r的值（可以用位运算加速），然后二分计算a^d mod n的值，最后把它平方r次。程序的代码比想像中的更简单，我写一份放在下边。虽然我已经转C了，但我相信还有很多人看不懂C语言。我再写一次Pascal吧。函数IsPrime返回对于特定的底数a，n是否是能通过测试。如果函数返回False，那说明n不是素数；如果函数返回True，那么n极有可能是素数。注意这个代码的数据范围限制在longint，你很可能需要把它们改成int64或高精度计算。function pow( a, d, n:longint ):longint;
begin
if d=0 then exit(1)
else if d=1 then exit(a)
else if d and 1=0 then exit( pow( a*a mod n, d div 2, n) mod n)
else exit( (pow( a*a mod n, d div 2, n) * a) mod n);
end;
function IsPrime( a,n:longint ):boolean;
var
d,t:longint;
begin
if n=2 then exit(true);
if (n=1) or (n and 1=0) then exit(false);
d:=n-1;
while d and 1=0 do d:=d shr 1;
t:=pow( a, d, n );
while ( d<>n-1 ) and ( t<>1 ) and ( t<>n-1 ) do
begin
t:=(t * t)mod n;
d:=d shl 1;
end;
exit( (t=n-1) or (d and 1=1) );
end;
对于大数的素性判断，目前Miller-Rabin算法应用最广泛。一般底数仍然是随机选取，但当待测数不太大时，选择测试底数就有一些技巧了。比如，如果被测数小于4 759 123 141，那么只需要测试三个底数2, 7和61就足够了。当然，你测试的越多，正确的范围肯定也越大。如果你每次都用前7个素数(2, 3, 5, 7, 11, 13和17)进行测试，所有不超过341 550 071 728 320的数都是正确的。如果选用2, 3, 7, 61和24251作为底数，那么10^16内唯一的强伪素数为46 856 248 255 981。这样的一些结论使得Miller-Rabin算法在OI中非常实用。通常认为，Miller-Rabin素性测试的正确率可以令人接受，随机选取k个底数进行测试算法的失误率大概为4^(-k)。
Miller-Rabin算法是一个RP算法。RP是时间复杂度的一种，主要针对判定性问题。一个算法是RP算法表明它可以在多项式的时间里完成，对于答案为否定的情形能够准确做出判断，但同时它也有可能把对的判成错的（错误概率不能超过1/2）。RP算法是基于随机化的，因此多次运行该算法可以降低错误率。还有其它的素性测试算法也是概率型的，比如Solovay-Strassen算法。另外一些素性测试算法则需要预先知道一些辅助信息（比如n-1的质因子），或者需要待测数满足一些条件（比如待测数必须是2^n-1的形式）。前几年AKS算法轰动世界，它是第一个多项式的、确定的、无需其它条件的素性判断算法。当时一篇论文发表出来，题目就叫PRIMES is in P，然后整个世界都疯了，我们班有几个MM那天还来了初潮。算法主要基于下面的事实：n是一个素数当且仅当(x-a)^n≡(x^n-a) (mod n)。注意这个x是多项式中的未知数，等式两边各是一个多项式。举个例子来说，当a=1时命题等价于如下结论：当n是素数时，杨辉三角的第n+1行除两头的1以外其它的数都能被n整除。
------------------------------------------------------------------------------------------------------------------------------------------------
这儿的快速幂写丑了。。。
前面都没有问题，我们主要来研究exit( (t=n-1) or (d and 1=1) )这句话;
or后面的意思是 a ^d mod p=1 or p-1
前面的意思是 a ^d在反复平方的过程中某一次 mod p=p-1，为什么是p-1呢，为什么不再 or t=1呢？
事实上，如果 a ^d在反复平方的过程中某一次 mod p=1,那么x ^2 mod p=1?
那么x=1 or p-1，又因为x不等于1，因为我们在最开始检验了 a ^d mod p是否=1 or p-1，
如果不是的话，在反复平方的过程中一定会先出现p-1，而不是1，而且如果出现1的话，那么一定之前出现了p-1，否则一直往前推，a ^d mod p=1 这与进入while循环矛盾
所以，这样的算法是没有问题的
我的代码：
1 var p,n,m:int64; t:longint; 2 procedure init; 3 begin 4 readln(p); 5 end; 6 procedure mul(x,y:int64;var z:int64); 7 var tmp:int64; 8 begin 9 tmp:=0; 10 while y>0 do 11 begin 12 if odd(y) then tmp:=(tmp+x) mod p; 13 y:=y>>1; 14 x:=(x+x) mod p; 15 end; 16 z:=tmp; 17 end; 18 19 function check(x:int64):boolean; 20 var cs,y,z:int64; 21 j:longint; 22 begin 23 cs:=n;y:=1; 24 while cs>0 do 25 begin 26 if odd(cs) then mul(y,x,y); 27 cs:=cs>>1; 28 mul(x,x,x); 29 end; 30 z:=n; 31 while (z<>p-1) and (y<>1) and (y<>p-1) do 32 begin 33 mul(y,y,y); 34 z:=z<<1; 35 end; 36 exit((odd(z)) or (y=p-1)); 37 end; 38 39 function isprime(p:int64):boolean; 40 var i:longint; 41 begin 42 if (not(odd(p))) or (p=1) then exit(false); 43 n:=p-1;m:=0; 44 while n and 1=0 do n:=n>>1; 45 for i:=1 to 10 do 46 if not(check(random(p-2)+2)) then exit(false); 47 exit(true); 48 end; 49 procedure main; 50 begin 51 if isprime(p) then writeln(‘Yes‘) else writeln(‘No‘); 52 end; 53 begin 54 assign(input,‘input.txt‘);assign(output,‘output.txt‘); 55 reset(input);rewrite(output); 56 readln(t); 57 while t>0 do 58 begin 59 dec(t); 60 init; 61 main; 62 end; 63 close(input);close(output); 64 end. 65
标签：style blog http color os strong io 数据
原文地址：http://www.cnblogs.com/zyfzyf/p/3888980.html
