轉載自Matrix大牛的博客 把代碼翻譯成C++
http://www.matrix67.com/blog/archives/234
題目鏈接：
http://hihocoder.com/problemset/problem/1287
一個數是素數（也叫質數），當且僅當它的約數只有兩個——1和它本身。規定這兩個約數不能相同，因此1不是素數。對素數的研究屬於數論范疇，你可以 看到許多數學家沒事就想出一些符合某種性質的素數並稱它為某某某素數。整個數論幾乎就圍繞着整除和素數之類的詞轉過去轉過來。對於寫代碼的人來說，素數比 想像中的更重要，Google一下BigPrime或者big_prime你總會發現大堆大堆用到了素數常量的程序代碼。平時沒事時可以記一些素數下來以 備急用。我會選一些好記的素數，比如4567, 124567, 3214567, 23456789, 55566677, 1234567894987654321, 11111111111111111111111 (23個1)。我的手機號前10位是個素數。我的網站域名的ASCII碼連起來(77 97 116 114 105 120 54 55 46 99 111 109)也是個素數。
素數有很多神奇的性質。我寫5個在下面供大家欣賞。
1. 素數的個數無限多（不存在最大的素數）
證明：反證法，假設存在最大的素數P，那么我們可以構造一個新的數2 * 3 * 5 * 7 * … * P + 1（所有的素數乘起來加1）。顯然這個數不能被任一素數整除（所有素數除它都余1），這說明我們找到了一個更大的素數。
2. 存在任意長的一段連續數，其中的所有數都是合數（相鄰素數之間的間隔任意大）
證明：當0<a<=n時，n!+a能被a整除。長度為n-1的數列n!+2, n!+3, n!+4, …, n!+n中，所有的數都是合數。這個結論對所有大於1的整數n都成立，而n可以取到任意大。
3. 所有大於2的素數都可以唯一地表示成兩個平方數之差。
證明：大於2的素數都是奇數。假設這個 數是2n+1。由於(n+1)^2=n^2+2n+1，(n+1)^2和n^2就是我們要找的兩個平方數。下面證明這個方案是唯一的。如果素數p能表示成 a^2-b^2，則p=a^2-b^2=(a+b)(a-b)。由於p是素數，那么只可能a+b=p且a-b=1，這給出了a和b的唯一解。
4. 當n為大於2的整數時，2^n+1和2^n-1兩個數中，如果其中一個數是素數，那么另一個數一定是合數。
證明：2^n不能被3整除。如果它被3除余1，那么2^n-1就能被3整除；如果被3除余2，那么2^n+1就能被3整除。總之，2^n+1和2^n-1中至少有一個是合數。
5. 如果p是素數，a是小於p的正整數，那么a^(p-1) mod p = 1。
這個證明就有點麻煩了。
首 先我們證明這樣一個結論：如果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小。
用同余式表述，我們證明了：
(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小定理(Fermat's Little Theorem)。Euler對這個定理進行了推廣，叫做Euler定理。Euler一生的定理太多了，為了和其它的“Euler定理”區別開來，有些地 方叫做Fermat小定理的Euler推廣。Euler定理中需要用一個函數f(m)，它表示小於m的正整數中有多少個數和m互素（兩個數只有公約數1稱 為互素）。為了方便，我們通常用記號φ(m)來表示這個函數（稱作Euler函數）。Euler指出，如果a和m互素，那么a^φ(m) ≡ 1 (mod m)。可以看到，當m為素數時，φ(m)就等於m-1（所有小於m的正整數都與m互素），因此它是Fermat小定理的推廣。定理的證明和Fermat小 定理幾乎相同，只是要考慮的式子變成了所有與m互素的數的乘積：m_1 * m_2 … m_φ(m) ≡ (a * m_1)(a * m_2) … (a * m_φ(m)) (mod m)。我為什么要順便說一下Euler定理呢？因為下面一句話可以增加我網站的PV：這個定理出現在了The Hundred Greatest Theorems里。
談到Fermat小定理，數學歷史上有很多誤解。很長一段時間里，人們都認為Fermat小定理的逆命題是正確的，並且有人親自驗證了 a=2, p<300的所有情況。國外甚至流傳着一種說法，認為中國在孔子時代就證明了這樣的定理：如果n整除2^(n-1)-1，則n就是素數。后來某個英 國學者進行考證后才發現那是因為他們翻譯中國古文時出了錯。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樹等方法使得查找偽素數表效率更高。
有 人自然會關心這樣一個問題：偽素數的個數到底有多少？換句話說，如果我只計算2^(n-1) mod n的值，事先不准備偽素數表，那么素性判斷出錯的概率有多少？研究這個問題是很有價值的，畢竟我們是OIer，不可能背一個長度上千的常量數組帶上考場。 統計表明，在前10億個自然數中共有50847534個素數，而滿足2^(n-1) mod n = 1的合數n有5597個。這樣算下來，算法出錯的可能性約為0.00011。這個概率太高了，如果想免去建立偽素數表的工作，我們需要改進素性判斷的算 法。
最簡單的想法就是，我們剛才只考慮了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，出錯的概率是否就可以降到0呢？沒想 到的是，居然就有這樣的合數，它可以通過所有a的測試（這個說法不准確，詳見我在地核樓層的回復）。Carmichael第一個發現這樣極端的偽素數，他 把它們稱作Carmichael數。你一定會以為這樣的數一定很大。錯。第一個Carmichael數小得驚人，僅僅是一個三位數，561。前10億個自 然數中Carmichael數也有600個之多。Carmichael數的存在說明，我們還需要繼續加強素性判斷的算法。
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或高精度計算（java）。
1 #include <iostream>
2 using namespace std ; 3 #define rd(x) (rand()%(x))
4 typedef unsigned long long ll; 5
6 ll pow_mod(ll a,ll b,ll r) 7 { 8 ll ans=1,buff=a; 9 while(b) 10 { 11 if(b&1) 12 ans=(ans*buff)%r; 13 buff=(buff*buff)%r; 14 b>>=1; 15 } 16 return ans; 17 } 18
19 bool test(ll n,ll a,ll d) 20 { 21 if(n==2) return true; 22 if(n==a) return false; 23 if(!(n&1)) return false; 24 while(!(d&1)) d>>=1; 25 ll t = pow_mod(a,d,n); 26 while(d!=n-1&&t!=n-1&&t!=1){ 27 t = t*t%n;//下面介紹防止溢出的辦法，對應數據量為10^18次方； 28 d<<=1; 29 } 30 return t == n-1||(d&1)==1;//要么t能變成n-1，要么一開始t就等於1 31 } 32
33 bool isprime(ll n) 34 { 35 int a[] = {2,3,5,7};//或者自己生成[2,N-1]以內的隨機數rand（）%（n-2）+2 36 for(int i = 0; i <= 3; ++i){
39 if(n==a[i]) return true; 40 if(!test(n,a[i],n-1)) return false; 41 } 42 return true; 43 } 44 int main() 45 { 46 int t,ans=0; 47 ll N; 48 for(cin >> t;t;t--){ 49 cin >> N; 50 cout << ((isprime(N))?"Yes":"No") <<endl; 51 } 52 }
對 於大數的素性判斷，目前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)。
后記：
注意上述算法中的冪運算是longlong類型，longlong×longlong肯定會出現溢出現象，如果不會java大整數，手里也沒有大整數乘法的模板
有一個小技巧可以避免溢出，方法就是乘法改為加法，把上面的代碼：
ll pow_mod(ll a,ll b,ll r) 7 { 8 ll ans=1,buff=a; 9 while(b) 10 { 11 if(b&1) 12 ans=(ans*buff)%r; 13 buff=(buff*buff)%r; 14 b>>=1; 15 } 16 return ans; 17 }
改為：
long long mod_mul(long long a, long long b, long long n) { long long res = 0; while (b) { if(b & 1) res = (res + a) % n; a = (a + a) % n; b >>= 1; } return res; } long long pow_mod(long long a, long long b, long long n) { long long res = 1; while(b) { if(b & 1) res = mod_mul(res, a, n); a = mod_mul(a, a, n); b >>= 1; } return res; }
就可以避免大整數，當然，復雜度也相應提高了一些。