定義 過去 lean 現在 我們 targe tro 不同的 greatest

以前判斷一個數n是否是素數我們一般都是暴力枚舉，O(sqrt(n))的算法。

最近看見了Miller Rabin算法，能以對數復雜度檢驗一個數是否為素數。

這個算法基於費馬小定理和二次探測定理。

下面轉載一個大佬 Matrix67 的博客：（原博客傳送門）

（華麗的分割線~）

---

數論部分第一節：素數與素性測試

一個數是素數（也叫質數），當且僅當它的約數只有兩個——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)也是個素數。還有，我的某個MM的八位生日也是一個素數。每次寫Hash函數之類的東西需要一個BigPrime常量時我就取她的生日，希望她能給我帶來好運。偶爾我叫她素MM，沒人知道是啥意思，她自己也不知道。
素數有很多神奇的性質。我寫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或高精度計算。
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整除。

Matrix67原創
轉貼請註明出處

---

（華麗的分割線No.2~）

再貼一下自己的代碼：

洛谷P3383 【模板】線性篩素數

在具體應用的過程中可以根據數據範圍選擇用來測試的素底數。上面大佬的文章也提到了。

 1 #include<cstdio>
 2 #define ll long long
 3 
 4 ll ksm(ll b,ll p,ll mod)
 5 {
 6     ll ret=1;
 7     while(p)
 8     {
 9         if(p&1)ret=(ret*b)%mod;
10         b=(b*b)%mod;
11         p>>=1;
12     }
13     return ret;
14 }
15 
16 ll pr[]={2,3,5,7,11,13,17,19};
17 
18 bool mrpt(ll n)
19 {
20     if (n==2)return true;
21     if((n&1)==0||n==1)return false;
22     for(ll i=0;i<8;i++)if(n==pr[i])return true;
23     ll temp=n-1,t=0,nxt;
24     while((temp&1)==0)temp>>=1,t++;
25     for(ll i=0;i<8;i++)
26     {
27         ll now=ksm(pr[i],temp,n);
28         nxt=now;
29         for(ll j=1;j<=t;j++)
30         {
31             nxt=(now*now)%n;
32             if(nxt==1&&now!=n-1&&now!=1)return false;
33             now=nxt;
34         }
35         if(now!=1)return false;
36     }
37     return true;
38 }
39 
40 ll n,m,opt;
41 
42 int main()
43 {
44     scanf("%lld%lld",&n,&m);
45     while(m--)
46     {
47         scanf("%lld",&opt);
48         if(mrpt(opt))printf("Yes\n");
49         else printf("No\n");
50     }
51     return 0;
52 }

Miller-Rabin素性測試