1. 根據(jù)概念判斷:

如果一個(gè)正整數(shù)只有兩個(gè)因子, 1和p，則稱p為素?cái)?shù).這是通常最先想到的判斷方法：
代碼:

bool isPrime(int n) {
if(n < 2) return false;
for(int i = 2; i < n; ++i)
if(n%i == 0) return false; return true; }

時(shí)間復(fù)雜度O(n).

2. 改進(jìn), 偶數(shù)都可以被2整除，去掉所有偶數(shù)

代碼：

bool isPrime(int n)
{
if(n < 2) return false;
if(n == 2) return true;
if(n%2 == 0) return false;

for(int i = 3; i < n; i += 2)
if(n%i == 0) return false;
return true;
}

時(shí)間復(fù)雜度O(n/2), 速度提高一倍.

3. 進(jìn)一步減少判斷的范圍

定理: 如果n不是素?cái)?shù), 則n有滿足1<d<=sqrt(n)的一個(gè)因子d.
證明: 如果n不是素?cái)?shù), 則由定義n有一個(gè)因子d滿足1<d<n.
如果d大于sqrt(n), 則n/d是滿足1<n/d<=sqrt(n)的一個(gè)因子.

代碼:

bool isPrime(int n)
{
if(n < 2) return false;
if(n == 2) return true;
if(n%2 == 0) return false;
int foo = (int)sqrt(n);
for(int i = 3; i <= foo; i += 2)
if(n%i == 0) return false;
return true;
}

時(shí)間復(fù)雜度O(sqrt(n)/2), 速度提高O((n-sqrt(n))/2).

4. 剔除因子中的重復(fù)判斷.
例如: 11%3 != 0 可以確定 11%(3*i) != 0.

定理: 如果n不是素?cái)?shù), 則n有滿足1<d<=sqrt(n)的一個(gè)"素?cái)?shù)"因子d.
證明: I1. 如果n不是素?cái)?shù), 則n有滿足1<d<=sqrt(n)的一個(gè)因子d.
I2. 如果d是素?cái)?shù), 則定理得證, 算法終止.
I3. 令n=d, 并轉(zhuǎn)到步驟I1.

由于不可能無(wú)限分解n的因子, 因此上述證明的算法最終會(huì)停止.

代碼: // primes[i]是遞增的素?cái)?shù)序列: 2, 3, 5, 7, ...
// 更準(zhǔn)確地說(shuō)primes[i]序列包含1->sqrt(n)范圍內(nèi)的所有素?cái)?shù)
bool isPrime(int primes[], int n)
{
if(n < 2) return false;
for(int i = 0; primes[i]*primes[i] <= n; ++i)
if(n%primes[i] == 0) return false;
return true;
} 假設(shè)n范圍內(nèi)的素?cái)?shù)個(gè)數(shù)為PI(n), 則時(shí)間復(fù)雜度O(PI(sqrt(n))).

函數(shù)PI(x)滿足素?cái)?shù)定理: ln(x)-3/2 < x/PI(x) < ln(x)-1/2, 當(dāng)x >= 67時(shí).

因此O(PI(sqrt(n)))可以表示為O(sqrt(x)/(ln(sqrt(x))-3/2)),

O(sqrt(x)/(ln(sqrt(x))-3/2))也是這個(gè)算法的空間復(fù)雜度.
5. 構(gòu)造素?cái)?shù)序列primes[i]: 2, 3, 5, 7, ...

由4的算法我們知道, 在素?cái)?shù)序列已經(jīng)被構(gòu)造的情況下, 判斷n是否為素?cái)?shù)效率很高;

但是, 在構(gòu)造素?cái)?shù)序列本身的時(shí)候, 是否也可是達(dá)到最好的效率呢?

事實(shí)上這是可以的! -- 我們?cè)跇?gòu)造的時(shí)候完全可以利用已經(jīng)被構(gòu)造的素?cái)?shù)序列!

假設(shè)我們已經(jīng)我素?cái)?shù)序列: p1, p2, .. pn

現(xiàn)在要判斷pn+1是否是素?cái)?shù), 則需要(1, sqrt(pn+1)]范圍內(nèi)的所有素?cái)?shù)序列,

而這個(gè)素?cái)?shù)序列顯然已經(jīng)作為p1, p2, .. pn的一個(gè)子集被包含了!

代碼: // 構(gòu)造素?cái)?shù)序列primes[]
void makePrimes(int primes[], int num)
{
int i, j, cnt;
primes[0] = 2;
primes[1] = 3;
for(i = 5, cnt = 2; cnt < num; i += 2)
{
int flag = true;
for(j = 1; primes[j]*primes[j] <= i; ++j)
{
if(i%primes[j] == 0)
{
flag = false;
break;
}
}
if(flag)
primes[cnt++] = i;
}
}
makePrimes的時(shí)間復(fù)雜度比較復(fù)雜, 而且它只有在初始化的時(shí)候才被調(diào)用一次.

在一定的應(yīng)用范圍內(nèi), 我們可以把近似認(rèn)為makePrimes需要常數(shù)時(shí)間.

在后面的討論中, 我們將探討一種對(duì)計(jì)算機(jī)而言更好的makePrimes方法.

6. 更好地利用計(jì)算機(jī)資源...

當(dāng)前的主流PC中, 一個(gè)整數(shù)的大小為2^32. 如果需要判斷2^32大小的數(shù)是否為素?cái)?shù),

則可能需要測(cè)試[2, 2^16]范圍內(nèi)的所有素?cái)?shù)(2^16 == sqrt(2^32)).

由4中提到的素?cái)?shù)定理我們可以大概確定[2, 2^16]范圍內(nèi)的素?cái)?shù)個(gè)數(shù).

由于2^16/(ln(2^16)-1/2) = 6138, 2^16/(ln(2^16)-3/2) = 6834,

我們可以大概估計(jì)出[2, 2^16]范圍內(nèi)的素?cái)?shù)個(gè)數(shù)6138 < PI(2^16) < 6834.

在對(duì)[2, 2^16]范圍內(nèi)的素?cái)?shù)進(jìn)行統(tǒng)計(jì), 發(fā)現(xiàn)只有6542個(gè)素?cái)?shù):

p_6542: 65521, 65521^2 = 4293001441 < 2^32, (2^32 = 4294967296)
p_6543: 65537, 65537^2 = 4295098369 > 2^32, (2^32 = 4294967296)

在實(shí)際運(yùn)算時(shí)unsigned long x = 4295098369;將發(fā)生溢出, 為131073.

在程序中, 我是采用double類型計(jì)算得到的結(jié)果.

分析到這里我們可以看到, 我們只需要緩沖6543個(gè)素?cái)?shù), 我們就可以采用4中的算法

高效率地判斷[2, 2^32]如此龐大范圍內(nèi)的素?cái)?shù)!

(原本的2^32大小的問題規(guī)模現(xiàn)在已經(jīng)被減小到6543規(guī)模了!)

雖然用現(xiàn)在的計(jì)算機(jī)處理[2, 2^16]范圍內(nèi)的6542個(gè)素?cái)?shù)已經(jīng)沒有一點(diǎn)問題,

雖然makePrimes只要被運(yùn)行一次就可以, 但是我們還是考慮一下是否被改進(jìn)的可能?!

我想學(xué)過java的人肯定想把makePrimes作為一個(gè)靜態(tài)的初始化實(shí)現(xiàn), 在C++中也可以

模擬java中靜態(tài)的初始化的類似實(shí)現(xiàn):

#define NELEMS(x) ((sizeof(x)) / (sizeof((x)[0])))

static int primes[6542+1];
static struct _Init { _Init(){makePrimes(primes, NELEMS(primes);} } _init;

如此, 就可以在程序啟動(dòng)的時(shí)候自動(dòng)掉用makePrimes初始化素?cái)?shù)序列.

但, 我現(xiàn)在的想法是: 為什么我們不能在編譯的時(shí)候調(diào)用makePrimes函數(shù)呢?

完全可以!!! 代碼如下:

代碼: // 這段代碼可以由程序直接生成
const static int primes[] = {
2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,
107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,
223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,
337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,
457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,
593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,
719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,
857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991, ... 65521, 65537 }; 有點(diǎn)不可思議吧, 原本makePrimes需要花費(fèi)的時(shí)間復(fù)雜度現(xiàn)在真的變成O(1)了!

7. 二分法查找

現(xiàn)在我們緩存了前大約sqrt(2^32)/(ln(sqrt(2^32)-3/2))個(gè)素?cái)?shù)列表, 在判斷2^32級(jí)別的

素?cái)?shù)時(shí)最多也只需要PI(sqrt(2^32))次判斷(準(zhǔn)確值是6543次), 但是否還有其他的方式判斷呢?

當(dāng)素?cái)?shù)比較小的時(shí)候(不大于2^16), 是否可以直接從緩存的素?cái)?shù)列表中直接查詢得到呢?

答案是肯定的! 由于primes是一個(gè)有序的數(shù)列, 因此我們當(dāng)素?cái)?shù)小于2^16時(shí), 我們可以直接

采用二分法從primes中查詢得到(如果查詢失敗則不是素?cái)?shù)).

代碼: // 缺少的代碼請(qǐng)參考前邊
#include <stdlib.h>
static bool cmp(const int *p, const int *q)
{
return (*p) - (*q);
}
bool isPrime(int n)
{
if(n < 2) return false;
if(n == 2) return true;
if(n%2 == 0) return false;

if(n >= 67 && n <= primes[NELEMS(primes)-1])
{
return NULL != bsearch(&n, primes, NELEMS(primes), sizeof(n), cmp);
}
else
{
for(int i = 1; primes[i]*primes[i] <= n; ++i)
if(n%primes[i] == 0) return false;
return true;
}
}
時(shí)間復(fù)雜度:

if(n <= primes[NELEMS(primes)-1] && n >= 67): O(log2(NELEMS(primes))) < 13;
if(n > primes[NELEMS(primes)-1]): O(PI(sqrt(n))) <= NELEMS(primes).

8. 素?cái)?shù)定理+2分法查找

在7中, 我們對(duì)小等于primes[NELEMS(primes)-1]的數(shù)采用2分法查找進(jìn)行判斷.

我們之前針對(duì)2^32緩沖的6453個(gè)素?cái)?shù)需要判斷的次數(shù)為13次(log2(1024*8) == 13).

對(duì)于小的素?cái)?shù)而言(其實(shí)就是2^16范圍只內(nèi)的數(shù)), 13次的比較已經(jīng)完全可以接受了.

不過根據(jù)素?cái)?shù)定理: ln(x)-3/2 < x/PI(x) < ln(x)-1/2, 當(dāng)x >= 67時(shí), 我們依然

可以進(jìn)一步縮小小于2^32情況的查找范圍(現(xiàn)在是0到NELEMS(primes)-1范圍查找).

我們需要解決問題是(n <= primes[NELEMS(primes)-1):

如果n為素?cái)?shù), 那么它在素?cái)?shù)序列可能出現(xiàn)的范圍在哪?

---- (n/(ln(n)-1/2), n/(ln(n)-3/2)), 即素?cái)?shù)定理!

上面的代碼修改如下:

代碼: bool isPrime(int n)
{
if(n < 2) return false;
if(n == 2) return true;
if(n%2 == 0) return false;

int hi = (int)ceil(n/(ln(n)-3/2));
if(n >= 67 && hi < NELEMS(primes))
{
int lo = (int)floor(n/(ln(n)-1/2));
return NULL != bsearch(&n, primes+lo, hi-lo, sizeof(n), cmp);
}
else
{
for(int i = 1; primes[i]*primes[i] <= n; ++i)
if(n%primes[i] == 0) return false;
return true;
}
} 時(shí)間復(fù)雜度:

if(n <= primes[NELEMS(primes)-1] && n >= 67): O(log2(hi-lo))) < ???;
if(n > primes[NELEMS(primes)-1]): O(PI(sqrt(n))) <= NELEMS(primes).


9. 打包成素?cái)?shù)庫(kù)(給出全部的代碼)

到目前為止, 我已經(jīng)給出了我所知道所有改進(jìn)的方法(如果有人有更好的算法感謝告訴我).

這里需要強(qiáng)調(diào)的一點(diǎn)是, 這里討論的素?cái)?shù)求法是針對(duì)0-2^32范圍的數(shù)而言, 至于像尋找

成百上千位大小的數(shù)不在此討論范圍, 那應(yīng)該算是純數(shù)學(xué)的內(nèi)容了.

代碼保存在2個(gè)文件: prime.h, prime.cpp. 代碼: // file: prime.h
#ifndef PRIME_H_2006_10_27_
#define PRIME_H_2006_10_27_
extern int Prime_max(void); // 素?cái)?shù)序列的大小
extern int Prime_get (int i); // 返回第i個(gè)素?cái)?shù), 0 <= i < Prime_max
extern bool Prime_test(int n); // 測(cè)試是否是素?cái)?shù), 1 <= n < INT_MAX
#endif
///
// file: prime.cpp
#include <assert.h>
#include <limits.h>
#include <math.h>
#include <stdlib.h>
#include "prime.h"
// 計(jì)算數(shù)組的元素個(gè)數(shù)
#define NELEMS(x) ((sizeof(x)) / (sizeof((x)[0]))) // 素?cái)?shù)序列, 至少保存前6543個(gè)素?cái)?shù)!
static const int primes[] = {
2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103, 107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211, 223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331, 337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449, 457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587, 593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709, 719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853, 857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991, ... 65521, 65537 };

// bsearch的比較函數(shù)
static int cmp(const void *p, const void *q)
{
return (*(int*)p) - (*(int*)q);
}

// 緩沖的素?cái)?shù)個(gè)數(shù)
int Prime_max()
{
return NELEMS(primes);
}

// 返回第i個(gè)素?cái)?shù)
int Prime_get(int i)
{
assert(i >= 0 && i < NELEMS(primes));
return primes[i];
}

// 測(cè)試n是否是素?cái)?shù)
bool Prime_test(int n)
{
assert(n > 0);
if(n < 2) return false;
if(n == 2) return true;
if(!(n&1)) return false;
// 如果n為素?cái)?shù), 則在序列hi位置之前
int lo, hi = (int)ceil(n/(log(n)-3/2.0));
if(hi < NELEMS(primes))
{ // 確定2分法查找的范圍
// 只有n >= 67是才滿足素?cái)?shù)定理
if(n >= 67)
lo = (int)floor(n/(log(n)-1/2.0));
else
{
lo = 0; hi = 19;
}
// 查找成功則為素?cái)?shù)
return NULL != bsearch(&n, primes+lo, hi-lo, sizeof(n), cmp);
}
else
{ // 不在保存的素?cái)?shù)序列范圍之內(nèi)的情況
for(int i = 1; primes[i]*primes[i] <= n; ++i)
if(n%primes[i] == 0) return false;
return true;
}
}
10. 回顧, 以及推廣

到這里, 關(guān)于素?cái)?shù)的討論基本告一段落. 回顧我們之前的求解過程, 我們會(huì)發(fā)現(xiàn)

如果缺少數(shù)學(xué)的基本知識(shí)會(huì)很難設(shè)計(jì)好的算法; 但是如果一味地只考慮數(shù)學(xué)原理,

而忽律了計(jì)算機(jī)的本質(zhì)特征, 也會(huì)有同樣的問題.

一個(gè)很常見的例子就是求Fibonacci數(shù)列. 當(dāng)然方法很多, 但是在目前的計(jì)算機(jī)中

都沒有實(shí)現(xiàn)的必要!

因?yàn)镕ibonacci數(shù)列本身是指數(shù)增長(zhǎng)的, 32位的有符號(hào)整數(shù)所能表示的位置只有前46個(gè):

代碼: static const int Fibonacci[] = {
0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,
2584,4181,6765,10946,17711,28657,46368,75025,121393,196418,
317811,514229,832040,1346269,2178309,3524578,5702887,9227465,
14930352,24157817,39088169,63245986,102334155,165580141,267914296,
433494437,701408733,1134903170,1836311903,-1323752223 };
因此, 我只需要把前46個(gè)Fibonacci數(shù)保存到數(shù)組中就可以搞定了!

比如: F(int i) = {return Fibonacci[i];} 非常簡(jiǎn)單, 效率也非常好. 《新程序員》：云原生和全面數(shù)字化實(shí)踐50位技術(shù)專家共同創(chuàng)作，文字、視頻、音頻交互閱讀

總結(jié)

以上是生活随笔為你收集整理的素数生成算法的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。