小白如何編程判斷一個數是否是質數?

張益唐如何編程判斷一個數是否是質數?

是兩個問題……

原題主沒有說清楚主語,有幾個答案就變成取消小白來問hello world的了,但事實上有一部分真需要編程來判斷一個數是否是質數的……


判定素數的方法啊……答案取決於很重要的兩個條件:「需要判斷的數有多大」,以及「這個數有沒有特殊形式」。以下我會按照適用範圍從小到大來盤點一下現有的判定演算法。

註:本文的脈絡參考了The Prime Pages上的介紹: https://primes.utm.edu/prove/index.html

寫在前面

當問題的輸入是一個正整數 [公式] 的時候,一般認為輸入規模是 [公式] 這是因為我們實際的輸入是 [公式] 的二進位展開式。後文中提到「多項式演算法」的時候,指的就是運行時間是 [公式] 的多項式的演算法。

另外文中提到的絕大多數演算法的主要部分都是大量模 [公式] 的乘法;文中提到的時間複雜度都假設我們用原始的 [公式] 的方式做這個乘法。如果 [公式] 很大,那實際計算中就要用到更優化的乘法方式(Karatsuba,Toom-Cook,包括FFT),整個演算法的時間複雜度也要相應的作調整。

naive的試除法

適用範圍: [公式] 很小,比如32位整數變數什麼的

對特殊形式的要求:無

時間複雜度: [公式]

普通的試除應該不用多加介紹了,從2試到 [公式] 就可以……可以用Eratosthenes篩法預先構造一個 [公式] 以內的素數表,用來優化試除過程。這個過程的缺點在於 [公式] 只要稍微大一點的時候,運行時間就是個天文數字了。

Fermat小定理

適用範圍:一切正整數

限制:由於偽素數的存在,Fermat小定理本身無法「證明」一個整數是素數

時間複雜度: [公式]

演算法本身大家也應該很熟悉了:挑一個底數 [公式] ,去計算 [公式] 是否成立。在使用快速冪演算法的情況下,這個計算過程只需要 [公式] 次模 [公式] 的乘法,而每次乘法的開銷是 [公式]

現在的問題在於,有的合數 [公式] 也會滿足Fermat小定理(這種n一般叫做「偽素數」),比如說 [公式] 更有甚者,有一類叫做Carmichael數的合數(最初的幾個例子是561,1105,1729,……),對所有的底數 [公式] (只要 [公式] 不是 [公式] 的因子)都滿足Fermat小定理。這就意味著Fermat小定理本身是無法證明 [公式] 是素數的

儘管如此,如果一個正整數 [公式] 通過了Fermat小定理的測試,那它是素數的可能性就很大了: [公式] 以內有455052511個素數,但是隻有14884個以2為底的偽素數,1547個Carmichael數。數學上一般把「通過了Fermat測試的正整數」稱為Probable Prime(PRP)。(這個術語好像還沒有通行的中文譯名)

一般面對一個較大的整數 [公式] ,在初步的試除過後,要做的下一件事就是用Fermat小定理(或者接下來要介紹的一些加強形式)來做出一個初步的判斷;如果 [公式] 通過了這個測試,說明它有極大的可能性是素數,然後可以用更高端的方法來證明 [公式] 確實是素數。

Fermat小定理的加強形式,Rabin-Miller測試,以及「工業級別」的素數

適用範圍:一切正整數

限制:仍然無法證明一個較大的整數是素數;但是隨機取 [公式] 個底數做測試的話,準確率可以達到至少 [公式]

時間複雜度: [公式]

Rabin-Miller測試是基於下面這個事實的:如果 [公式] 是素數,那麼 [公式] 的情況下,1的平方根只能是 [公式] 。把這個東西跟Fermat小定理結合起來,就能發現:

  • [公式]
  • 如果上面一個式子的右端是1,那麼就可以繼續開平方根得到 [公式]
  • 依次類推,直到左邊的指數 [公式] 是奇數不能再開平方根,或者某一步右邊變成了-1為止。
  • 如果某一步開平方根的時候右端直接從1變成了「不是 [公式] 的數」,則 [公式] 是合數。

這個方法雖然也有偽素數的存在(比如 [公式] 的時候2047可以通過Rabin-Miller測試),但是不會有類似Carmichael數的東西出現了:對於任意的合數 [公式] , 總存在一個適當的底數 [公式] 使得 [公式] 無法通過以 [公式] 為底的Rabin-Miller測試。事實上,在 [公式] 這個範圍裏,至少有75%的底數滿足這個條件。這就意味著,如果我們以隨機的方式取出 [公式] 個底數 [公式] ,而且正整數 [公式] 通過了所有的以 [公式] 為底的Rabin-Miller測試,那麼從某種意義上說, [公式] 不是素數的概率最多是 [公式] 。舉個例子,如果[公式] ,那麼這個概率就是 [公式] ,幾乎可以忽略了。所以Rabin-Miller測試是在實際的密碼學應用中最常用的判定方法,通過了這個測試的數可以看做是「工業級別」的素數。

額外說一句:很多數學軟體裏的isprime()或者類似的函數,使用了底分別是2,3,5,7,11,13,17的7次Rabin-Miller測試。這個測試當然不適用於所有的正整數,但是它的最小反例是 341550071728321……所以對於這個值以下的整數來說,「通過這7次Rabin-Miller測試」就可以證明它是素數了。

再額外說一句:如果我們假設廣義Riemann猜想(GRH)成立,那麼證明 [公式] 是素數就只需要驗證「 [公式] 可以通過以 [公式] 為底的Rabin-Miller測試」了。這個演算法,如果能嚴格證明的話,運行時間是 [公式]

Lucas序列

適用範圍:一切正整數

限制:和Fermat小定理一樣,無法「證明」一個整數是素數

時間複雜度:也是 [公式]

大家可能聽說過Fibonacci數列滿足這麼一條性質:如果 [公式] 是素數,則 [公式] ,其中 [公式] 是二次剩餘的Jacobi符號。這其實是二次擴域 [公式] 裡面的Fermat小定理的推廣:

[公式]

這個結果自然也可以推廣到別的二次擴域上:

[公式]

實際計算的時候就會涉及Fibonacci序列的類似物:一般的有齊次線性二階遞推關係的序列(叫做Lucas Sequence)。這種序列一般的定義是

[公式]

也有通項公式

[公式] ,其中 [公式]

這時需要驗證的同餘式就變成了 [公式]

[公式] 也可以用類似於快速冪的方式進行計算,所以驗證這個同餘式的時間複雜度也是[公式]

理所當然,和剛才提到的Fermat小定理一樣,這類測試也會有偽素數的存在(Fibonacci數的情況下,最小的例子是323和377),也有類似Rabin-Miller測試的強化版本。

一句題外話:把Rabin-Miller方法和強化的Lucas方法組合起來,我們就得到了當前最強大的Probable Prime判定方法:如果一個正整數 [公式] 通過了底是2的Rabin-Miller測試,也通過了底是 [公式] 的強化版Lucas測試(其中 [公式] 是序列 [公式] 之中第一個滿足 [公式] 的元素),那 [公式] 幾乎可以確定是素數。這一般叫做Baillie-PSW Test,雖然不嚴格但是至今為止還沒有發現反例!

-----------------------分割線-------------------

以上介紹了一些較快的可以「初步測試」一個正整數是不是素數的辦法。如果某個 [公式] 通過了上面的這些測試,那它肯定有很大的可能性是素數。接下來的問題是:如何去嚴格的證明 [公式] 就是素數呢?

接下來有請我們的各種高端方法出場(掌聲)。

n-1方法

適用範圍:有特殊形式的正整數 [公式]

對特殊形式的要求: [公式] 的素因子大部分已知

時間複雜度:[公式]

我們先來回顧一下對Fermat小定理的一個羣論證明:如果 [公式] 是素數,則乘法羣 [公式] 的階數是 [公式] ,由羣論中的Lagrange定理立得。

這個證明告訴了我們什麼?Fermat小定理包含了一些關於羣 [公式]的信息。如果我們能想辦法說明 [公式]的階數就是 [公式] ,那就足以證明 [公式] 是素數了。如果我們知道 [公式] 的所有素因子,那這個驗證過程並不困難。

定理(Lucas, 1891): 如果存在一個底數 [公式] ,使得以下兩條成立:

  • [公式]
  • [公式] 的每個素因子 [公式] ,都有 [公式]

那麼 [公式] 是素數。

證明:這兩個條件組合起來可以說明 [公式][公式] 中的階數就是 [公式] ,所以[公式]的階數恰好就是 [公式] ,從而 [公式] 是素數。

當然這個演算法本身也擁有很大的改進餘地。因為對於很多整數 [公式] 來說我們並不能保證知道 [公式] 的所有素因子(因子分解一般比素數判定難得多),我們接下來介紹一些「只知道 [公式] 的一部分因子」的情況下也能生效的辦法。

定理(Pocklington, 1914): 假設 [公式] ,其中 [公式] 的素因子都已知。如果存在一個底數 [公式] ,使得以下兩條成立:

  • [公式]
  • [公式] 的每個素因子 [公式] ,都有 [公式]

那麼 [公式] 的每個素因子都是 [公式] 的形式。特別地,如果 [公式] ,則 [公式] 是素數。

這個定理說明瞭把 [公式] 分解到「一半」就可以證明 [公式] 是素數了。後續的一些工作用更加細緻的分析把這個「一半」的界限降低到了三分之一 [Brillhart-Lehmer-Selfridge 1975], 30% [Konyagin-Pomerance 1997], 以及 [公式] [Coppersmith-Howgrave-Graham, 2008](後兩篇文章的演算法需要做一些附加的計算,涉及Lattice Reduce以及LLL演算法)。 篇幅所限,此處不做進一步展開,有興趣的讀者可以自行查閱文獻。

n+1方法

適用範圍:有特殊形式的正整數 [公式]

對特殊形式的要求: [公式] 的素因子大部分已知

時間複雜度:[公式]

[公式] 方法和 [公式] 方法非常類似,只不過涉及的羣變成了商羣 [公式] 。如果 [公式] 是素數且 [公式] ,那麼這個羣的階數就應該是 [公式] ;反過來,如果我們用跟上一節類似的方法驗證了這個羣的階數就是 [公式] ,那也就說明瞭 [公式] 是素數。這個驗證的過程用到了之前提到的Lucas序列:如果存在合適的 [公式] 使得 [公式] ,且相應的Lucas序列滿足

[公式] ,那麼 [公式] 就是素數。

上一節所說的改動同樣生效:如果 [公式][公式] 的素因子都已知且滿足

[公式] ,那麼 [公式] 的每個素因子都是 [公式] 的形式。同樣的,如果這個 [公式] 足夠大(和上一節的結果平行,最差只需要 [公式] ),那麼就可以通過一些附加的計算來證明 [公式] 是素數。

一句題外話:大家可能聽說過檢驗Mersenne素數 [公式] 的Lucas-Lehmer測試。這其實是 [公式] 方法的一個特殊情況:這裡 [公式] ,它的素因子自然只有2,而Lucas-Lehmer測試等價於底是 [公式] (或者說[公式] )的時候驗證了

[公式]

這兩節介紹的 [公式] 方法是1890-1980年這近一個世紀的時間裡僅有的可以證明某個很大的正整數確實是素數的方法。現在所知道的Top-5000的大素數(參考https://primes.utm.edu/primes/home.php)無一例外都是用這兩種方法之一證明的,它們的形式也都是「某個素因子都已知的數 [公式] 」這樣。

順帶一提:如果 [公式] 都有一些已知的素因子(比如 [公式] ),那分別驗證完兩部分條件之後只要 [公式] 就可以證明瞭。這一點對於「 [公式] 的素因子分解並不平凡,但是仍然可以分解出相當一部分」的這種 [公式] 最有效。典型的例子是Fibonacci素數: [公式][公式] ,可以由Fibonacci序列的整除性質得到 [公式] 的大量素因子。

----------------------70-80年代的分割線------------------------

以上的 [公式] 方法通常被稱為「古典的素數判定法」,共同要求是我們需要知道 [公式] 的一定數量的素因子。對於一般的沒有特殊形式的正整數 [公式] 來說,對 [公式] 做因子分解是一件很困難的事。以下介紹80年代以後開發出來的不需要 [公式] 的特殊形式的「現代素數判定法」。

APR-CL演算法

適用範圍:一切正整數

時間複雜度: [公式]

剛才介紹的 [公式] 方法的本質都在於「確定一個羣的大小」: [公式] 的時候是 [公式] ,而 [公式] 的時候涉及到 [公式] 的一個二次擴張。70年代末期,有人開始考慮更高次的擴張會不會對問題有幫助。這個方向最初的結果是Williams做出的:通過考慮次數是3,4,6的擴張,他們說明瞭 [公式] 以及 [公式] 的素因子也可以對證明 [公式] 是素數做出貢獻。然而問題來了:這些新增的素因子一般數量太少,對問題沒有本質的影響。

80年代的時候,有一批研究者決定一不做二不休,乾脆考慮更高次數的擴張:一般的 [公式] 次擴張可以讓我們利用 [公式] 的素因子(具體的計算會涉及分圓域裏的Gauss和)。這有什麼好處呢?好處在於:如果 [公式] 本身有很多素因子,那麼Fermat小定理保證了 [公式] 會包含很多「免費」的素因子:如果某個素數 [公式] 滿足 [公式] ,那就可以直接推出 [公式] ,這個過程已經和 [公式] 無關了。[Adleman-Pomerance-Rumely 1983]裡面提到了一個典型的例子:當 [公式] 的時候, [公式] 的「免費」素因子的乘積是

[公式]

這意味著對於 [公式] 以內的正整數來說,這些「免費」的 [公式] 的素因子的乘積就已經超過了 [公式] ,接下來就可以直接證明 [公式] 是素數了。(為何是 [公式] ?這是因為我們現在的計算是在一個很大的分圓域裏做的,這種情況下「之前提到的進一步的改進」會嚴重影響運行時間。) 如果將這個演算法一般化,有結果證明總可以找到一個 [公式] 使得相應的「免費素因子」的乘積超過 [公式] ,這就是開始提到的那個時間複雜度的來源。

Adleman,Pomerance 和 Rumely 提出這個演算法之後,[Cohen-Lenstra, 1984] 對其中的計算過程做出了重要的改進(粗略地說:用Jacobi和代替了Gauss和,這樣需要處理的分圓域的次數就大大降低了)。改進後的演算法就以這五個人的名字命名,叫做APR-CL演算法。這是歷史上第一個相對快速的可以對任意正整數生效的素數判定/證明演算法。

最後提一句:這個時間複雜度理論上不是多項式演算法,但是實際應用中 [公式] 的增長極慢,幾乎是個常數,所以實際的運行時間「幾乎是多項式的」。

橢圓曲線演算法(Elliptic Curve Primality Proving, ECPP)

適用範圍:一切正整數

時間複雜度: [公式] (有隨機性,這個是期望運行時間);演算法會生成相對簡短的證書(空間複雜度 [公式] ),通過這個證書可以快速驗證 [公式] 是素數(時間複雜度 [公式] )而不用重複整個計算過程

這個演算法可以看做是 [公式] 方法的另一方面的推廣:將乘法羣 [公式] 推廣到了一般的橢圓曲線羣 [公式] 。使用橢圓曲線羣的優點很明確:它的階數不再是固定的 [公式] ,而是一個隨著曲線變化,取值範圍可以落在 [公式] 這個區間裡面的整數。和上文中提到的Pocklinton定理類似,ECPP的理論基礎是如下的定理:

定理(可能已經是forklore了):如果存在一條 [公式] 上的橢圓曲線 [公式] ,以及曲線上的一個點 [公式] ,滿足以下幾個條件:

  • [公式] ,其中 [公式] 是素數且 [公式] ;
  • [公式] 中, [公式]

[公式] 是素數。

利用這個定理證明 [公式] 是素數,需要兩個先決條件:找到合適的橢圓曲線 [公式] ,以及證明我們用到的 [公式] 確實是素數。

第二個問題很好解決,因為我們把一個關於 [公式] 的問題化歸成了關於 [公式] 的問題,而一般來說 [公式] ,所以這個演算法可以無限遞歸下去,直到涉及的素數規模足夠小可以用Rabin-Miller甚至試除來搞定。剩下的就是第一個問題:怎樣找到合適的橢圓曲線,使得相應的羣的階數是 [公式] 的形式?

Goldwasser-Kilian的第一版ECPP演算法使用了非常原始的方式:隨機取一條橢圓曲線 [公式] (嚴格的說,隨機取 [公式] 中的 [公式][公式] ),用Schoof演算法算出階數 [公式] ,試著對算出來的階數做因子分解;如果得不到 [公式] 這個形式就扔掉換一條。演算法最大的硬傷在於Schoof演算法本身太慢( [公式] ),拖慢了整個ECPP方法的運行時間。

1993年的Atkin-Morain演算法解決了這個問題。他們注意到如果一條橢圓曲線有復乘(complex multiplication),那麼相應的曲線羣的階數就很好計算(這個理論背景我不太懂,希望有熟悉相關領域的朋友在評論區補充一下);所以只要把「隨機選擇一條曲線」的範圍限定在complex multiplication的範圍裏,就可以大大提高整個演算法的效率。

需要指出一點:整個計算過程中最費時間的步驟是「找到合適的橢圓曲線,並且對曲線羣的階數進行分解」。如果在運行過程中把這些信息記錄下來(準確的說:在這個過程中的每一步,記錄下來四元組 [公式] ),那它們就可以用來快速的驗證 [公式] 就是素數。這些記錄下來的信息就是本節開頭提到的證書(certificate)

ECPP演算法是目前為止最快速的「不需要特殊形式」的素數判定演算法。一般對於 [公式] 左右的素數,在當今普通的個人電腦上只需要不到半個小時就可以完成;當前(截止到2019年1月18日)用這個演算法證明的最大素數有 [公式] 個十進位位(參考https://primes.utm.edu/top20/page.php?id=27),據維基說運行時間大概是兩年。

AKS演算法

適用範圍:一切正整數

時間複雜度: [公式] ;空間複雜度: [公式] ;為啥在這要提到空間複雜度?因為上文提到的一切演算法的空間複雜度都是 [公式] ,只有AKS不一樣

注意:這個方法目前只有理論價值,沒有實用價值!

AKS演算法的本質可以看作「驗證特定的多項式環上的Fermat小定理」。它的大概過程如下:

  • 輸入一個正整數 [公式] 。檢驗 [公式] 不是其他正整數的冪。
  • 找到最小的正整數 [公式] 使得 [公式][公式] 中的階數至少是 [公式]
  • 一路試除到 [公式]
  • 對於所有的 [公式] ,驗證下面的多項式同餘式: [公式]
  • 如果都通過則 [公式] 是素數。

它的理論價值在於:這是第一個確定性的,對所有整數有效不依賴於其它猜想多項式時間演算法。跟上文中的幾個演算法作對比:

  • [公式] 方法需要用到 [公式] 的特殊形式
  • APR-CL演算法的時間複雜度不是多項式
  • Rabin-Miller方法依賴於GRH的成立
  • ECPP方法滿足上面提到的別的條件,但是涉及到「隨機選取橢圓曲線」的過程,所以不是確定性的。

然而………………………………這個理論上很美的演算法實際中並沒有什麼卵用。實機測試顯示,對於 [公式] 左右的正整數 [公式] ,APR-CL和ECPP都是秒出結果,而AKS的運行時間已經要以天計了。

最後總結一下

  • 如果要判斷的正整數 [公式] 很小(int32),直接試除就好。
  • 如果 [公式] 稍微大一點 [公式] ,用7次Rabin-Miller。
  • [公式] 再大一些的時候,首先用試除和Baillie-PSW做一個初步測試。如果通過了初步測試,再用以下的方法來證明 [公式] 是素數:
  • 如果 [公式] 很容易做因子分解,用相應的 [公式] 方法。
  • 如果 [公式] 沒有特殊形式,用APRCL或者ECPP。
  • 如果只是想滿足好奇心,可以試著實現一下AKS……

--------------------------------完結撒花--------------------------


又是喜聞樂見的素數判定問題。本蒟蒻來強答一發

一個人人皆知的方法

你只要知道質數的定義,就可以報出一個正確的 [公式] 做法;當然是從 [公式][公式] 一個一個枚舉,看看除了 [公式] 和它本身之外還有沒有其他質因數了!

當然這個方法顯然太慢了一點,想要快速判斷 [公式] 這種級別的質數怎麼辦?

改進方法有了!觀察一下它的質因數,如果是合數那麼必然存在一個 [公式] 的質因數(想想為什麼)。於是我們只要枚舉 [公式][公式] 內的所有數,看看是不是都不能整除,就行啦!

換一種角度

觀察一下費馬小定理。(這裡就不寫了,可以看費馬小定理_百度百科 )

發現如果它的逆命題成立,那事情就變得十分美好了---隨機一個值我們就可以知道 [公式] 是不是質數。然而現實總是殘酷的,費馬小定理的逆命題並不成立。但不妨礙我們設計出如下演算法:

[公式] 為奇數時(偶數時特判一下)讓 [公式][公式] 之間(包括 [公式][公式] )選取隨機值,如果等式不成立,那麼 [公式] 肯定不是素數,如果成立,那麼 [公式] 就有較大可能是素數,稱為偽素數

然而這種演算法誤判的概率實在是太高了(最典型的例子,Carmichael數),所以它並沒有什麼卵用。嘗試改進這種方法,於是就有了現在人們耳熟能詳的Miller-Rabin素數判定演算法,它基於二次探測定理(可看Miller_Rabin素數測試[Fermat小定理][二次探測定理][同餘式][Wilson定理] - pi9nc的專欄 - CSDN博客 ):

(版權糾紛不存在的啦,這是窩的課件)

可以發現這種演算法極大地提升了正確率。在 [公式] 範圍內出錯的只有一個數(忘記是什麼了,哪位大佬評論區說一下QwQ)

一般現在我們檢測素數都是用的Miller-Rabin啦~


八卦一下,其實Miller原先提出的是一種基於廣義黎曼猜想的 [公式] 演算法,然後被Rabin教授改成了不基於廣義黎曼猜想的隨機化演算法~(≧▽≦)/~

具體可以看

https://core.ac.uk/download/pdf/82649441.pdf?

core.ac.uk


Update 1.11 有答主提到了AKS演算法,目前的複雜度估計是 [公式] 的。假設Agrawals conjecture成立那麼它的複雜度將降至 [公式] 。然而測試速度還是遠遠不及Miller-Rabin(並且甚至更多的證據表明這個猜想並不成立),所以Miller-Rabin仍然是現在最廣泛採用的素數判定演算法。


update 1.17

昨天 @醬紫君 譴責我沒有講一些modern的做法

但可是我已經看過mma的幫助文檔了啊...

大家以後都去學ECPP,別來學米勒拉賓了(逃)


唔,當然了,醬紫早就學完代數數論了,所以嘛...


AKS演算法,可以在多項式時間內確定一個數是否為素數

PRIMES is in P | Annals of Mathematics?

annals.math.princeton.edu

原文給出的複雜度是 [公式] ,後來有其他文章改進到 [公式]


我見過的實現版本大致是這樣的

def is_prime(n):
不通過 = 用已經篩好的素數表試除(n)
if 不通過:
return False
# https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
不通過 = MillerRabin素數檢測(n)
if 不通過:
return False
不通過 = 檢測是否為非平方數(n)
if 不通過:
return False
# https://en.wikipedia.org/wiki/Lucas_pseudoprime
不通過 = StrongLucas素數檢測(n)
if 不通過:
return False
return True

這裡面的核心就是同時使用 Miller-Rabin 和 Lucas 來檢測,雖然兩者都有非常小的可能性產生漏判,但是由於兩者是獨立的,所以同時產生漏判的概率可以說是微乎其微,反正迄今為止應該還沒有找到一個反例。


matlab ,將就著看吧

a=input(請輸入一個整數(end in 0):);

while(a~=0)

if a==1

disp(1 is not a prime or a composite number.)

elseif isprime(a)==1

fprintf(%d is a prime.
,a)

elseif isprime(a)==0

f=factor(a);m=length(f);

fprintf(%d is not a prime,%d=%d,a,a,f(1))

for i=2:m

fprintf(*%d,f(i));

end

fprintf(
);

end

a=input(請輸入一個整數(end in 0):);

end


推薦閱讀:
相關文章