Fermat素性測試， Miller-Rabin素性測試
昨天看了看spoj的第2題（坑啊~）， 說的是找出 A,B之間的素數， 有T<=10組資料， A, B均小於10^9， A，B之差小於10^5。
打表是不可能了， 只能一個一個判斷。
這是我們需要一個強力的方法來判斷一個數是素數呢，還是合數。 所以可以想到這個素性測試。
Fermat素性測試
首先來說一個性質： 如果p是素數，a是小於p的正整數，那麼a^(p-1) mod p = 1。
這個定理的證明就不多講了（鍵盤的shift鍵快要敲爛了）， 請看第一句。
不過有一些合數也可以通過這樣的一個測試，當a=2時， 前10E有5597個合數。 所以我們可以試多一點數。（不過還是有反例， 請百度Carmichael數
若這個方法可以不讓合數漏網， 這個還是很容易實現的。
那麼對於一個數， 只要判斷一個素數a的p-1次方mod p的值是否為1。
對於每一個區間A,B， 先計算出a^A(a是一個數， a^A可以用快速冪)，再依次判斷即可。
時間複雜度 O（T * A *（log N + L）） A為需要列舉為a的數的個數（有點繞）。
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)。
所以我們可以得到一個判斷n是否為素數的方法：
原來判斷素數的式子為 a^(p-1) mod p = 1， 若 p-1 % 2==0 ， 那麼原式就變成了 ： (a ^ x) ^ 2 mod p = 1. (x = p/2)
若式子不成立， 那麼它連Fermat素性測試也沒有通過，則肯定不是素數， 那麼成立以後，
若p為素數， 根據上面的式子， a^x= 1 或 a^x= p-1.
所以有兩種情況： 當x=0【即再也分解不出2因子】或 a^x= p-1 時， 我們宣佈： 這個傢伙極有可能是素數。
所以實現方法如下：
不斷地提取指數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。
但遺憾的是，竟然還是有一些這樣的合數，能通過這樣的測試， 若a=2，通過測試的第一個合數是2407.
加上a=3， 第一個反例到了1 373 653。 這還是比較萬幸的,至少沒有能通過所有測試的合數
題目說數都是小於10E， 那麼用a=2, 3, 61是沒有問題的， 在4 759 123 141之前都是沒有問題的。
附上標程： 參見本文第一句：
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; // 提出因子2
t:=pow( a, d, n ); // 運用快速冪計算出a^d
while ( d<>n-1 ) and ( t<>1 ) and ( t<>n-1 ) do // 即還有因子2， mod 不等於 1 或 n-1
begin
t:=(t * t)mod n;
d:=d shl 1;
end;
exit( (t=n-1) or (d and 1=1) ); // 跳出 ， mod的數等於 n-1 或 一開始就成立。
end;所以這道題解決了。時間複雜度為O (N log N)。注意，時限為6s。