素?cái)?shù)的計(jì)算 - 從試除到篩法
昨天搜索素?cái)?shù)的問題的時(shí)候，找到一篇很棒的文章，轉(zhuǎn)載一下，并加上一些自己的理解。
文章鏈接: 素?cái)?shù)的計(jì)算: 從試除到篩法(C++實(shí)現(xiàn))
素?cái)?shù)定理：
素?cái)?shù)的個(gè)數(shù)是有規(guī)律的，對(duì)于正實(shí)數(shù)x，定義T(x)為素?cái)?shù)計(jì)數(shù)函數(shù)：不大于x的素?cái)?shù)的個(gè)數(shù)。
那么會(huì)有：T(x) 約等于 x / ln(x) .
其中 T(x) / (x/ln(x)) 總是小于1.17。這個(gè)性質(zhì)能夠幫我們確定x以內(nèi)的最大素?cái)?shù)個(gè)數(shù)，以分配存儲(chǔ)空間(即數(shù)組開多大)。
對(duì)于找不大于x的所有素?cái)?shù)，有很多方法，下面從最簡(jiǎn)單的到最高效的介紹。
試除法
思想：
判斷方法：對(duì)于一個(gè)數(shù)n，將其分別除以[2,n]的每一個(gè)整數(shù)，如果有可以整除的，則不是素?cái)?shù)。
循環(huán)所有不大于x的數(shù)，確定數(shù)xi是否為素?cái)?shù)以求出不大于x的所有素?cái)?shù)。
算法實(shí)現(xiàn):
#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
using namespace std;
// 判斷n是否為素?cái)?shù)
bool is_prime(int n)
{
if (n < 2)
return false;
for (int i = 2; i < n; ++i)
if (n % i == 0)
return false;
return true;
}
// 計(jì)算所有不大于n的素?cái)?shù)
void get_prime(vector<int>& prime, int n)
{
for(int i = 2; i <= n; ++i)
if(is_prime(i)) // 判斷i是否是素?cái)?shù)
prime.push_back(i);
}
int main()
{
int n = 100000;
vector<int> prime;
get_prime(prime, n);
return 0;
}
算法復(fù)雜度: O(n)
試除法優(yōu)化
思想：
對(duì)上面思想的分析，我們發(fā)現(xiàn)，若有數(shù)n = x*y，那么y與x中必有一個(gè)滿足：k <= sqrt(n)。即在sqrt(n)之前沒有找到能夠整數(shù)n的數(shù)，那么n之后也不會(huì)有，所以對(duì)于判斷n是否為素?cái)?shù)的方法中的循環(huán)，只需要判斷到sqrt(n)即可。
算法實(shí)現(xiàn)
更改is_prime方法：
bool is_prime(int n)
{
if (n < 2)
return false;
for (int i = 2; i * i <= n; ++i)
if (n % i == 0)
return false;
return true;
}
算法復(fù)雜度： O(n * sqrt(n))
上面的都是暴力手法，下面介紹科學(xué)的藝術(shù)手法。
埃氏篩法
埃拉托斯特尼篩法，是古希臘數(shù)學(xué)家發(fā)明的計(jì)算素?cái)?shù)的方法(媽耶，古希臘就研究出來了)。
思想：
對(duì)于求解不大于n的素?cái)?shù)：
找出不大于sqrt(n)內(nèi)的素?cái)?shù) p1、p2、p3 .... pn (1 <= n <= sqrt(k))
依次剔除不大于n的pi的倍數(shù)。
剩下的都是素?cái)?shù)。
這里有一張埃氏篩法的工作原理圖：
算法實(shí)現(xiàn)：
void get_prime (vector<bool> &isPrime, int n) {
isPrime.assign(n + 1, true);
if (n < 2) return;
for (int i = 2; i <= n; i++) {
if (isPrime[i]) {
// 若計(jì)算規(guī)模較大，加上 i < sqrt(n) 的判斷條件
for (int j = i * i; j <= n && i < sqrt(n); j += i) {
isPrime[j] = false;
}
}
}
}
這里有一個(gè)原文沒有提到的點(diǎn)，在我的實(shí)踐中，當(dāng) n > 46000 時(shí)，會(huì)出現(xiàn) i*i 超出int范圍的問題，此時(shí)j會(huì)是負(fù)值，造成非法數(shù)組訪問的錯(cuò)誤。所以如果計(jì)算規(guī)模較大，加上 i < sqrt(n) 的判斷條件。
算法時(shí)間復(fù)雜度： O(nloglogn)
歐拉篩法
思想
埃氏篩法會(huì)對(duì)兩個(gè)素?cái)?shù)的公倍數(shù)多次剔除，根據(jù)這一問題優(yōu)化后，便出現(xiàn)了歐拉篩法：對(duì)于一個(gè)沒有被篩過的數(shù)，只要被第一個(gè)數(shù)篩了就行了。
算法實(shí)現(xiàn)
void get_prime(vector<int> &prime, vector<bool> &isPrime, int n) {
isPrime.assign(n + 1, true);
int max_num = (n / log(n)) * 1.17 + 1;
if (n < 2) return;
for (int i = 2; i <= n; i++) {
if (isPrime[i]) {
prime.push_back(i);
}
for (int j = 0; j < isPrime.size() && i * prime[j] <= n; j++) {
isPrime[i*prime[j]] = false;
if (i % prime[j] == 0) break;
}
}
}
這個(gè)算法的精髓主要在于：
if (i % prime[j] == 0) break;
證明：對(duì)于一個(gè)數(shù)i，若能整除prime[j]，那么有 i = a * prime[j]。當(dāng)判斷prime[j+1]能否被i整除時(shí)，有 i * prime[j + 1] = a * prime[j] * prime[j+1]，即此數(shù)(i*prime[j+1])已經(jīng)被prime[j]剔除過了，所以跳出。
此算法時(shí)間復(fù)雜度：O(n)
數(shù)組開多大
最前面已經(jīng)說了，數(shù)組開多大可以用T(x) / (x/ln(x)) <= 1.17 來決定，那么：
int max_num = (n / log(n)) * 1.17 + 1;
