篩素數，判素數，1E11以内素數個數
篩素數
埃式篩法O(nlog(log(n)) )
const int N = 1e7+2;
bool notPrime[N];
int sushu[700000],cnt;
void Ero_getPrime() { //埃式篩法
int n = N-1;
cnt = 0;
for(int i=2;i<=n;++i) {
if(!notPrime[i]) {
sushu[cnt++] = i;
if(i<=n/i)
for(int j=i*i;j<=n;j+=i) //i*i之前的合數都被更小的質因數篩過了
notPrime[j] = true;
}
}
}
2篩掉了4, 6, 8, 10, 12,14,16,18,20,22,24,26,28,30,32…
3篩掉了9,12,15,18,21,24,27,30,33…
5篩掉了25,30,35,50,45…
我們發現：30 = 2*15 = 3*10 = 5*6, 被篩了不止一次。
j從i*i開始，是因為i*i之前的合數，一定被它最小的質因子篩過了(可能也被更大的質因子篩過)
歐拉篩(線性篩)。每個合數被其最小的質因子篩掉，隻篩一次，接近O(n)
const int N = 1e7+2;
bool notPrime[N];
int sushu[700000],cnt;
void getPrime() {
cnt = 0;
int n = N-1;
for(int i=2;i<=n;++i) {
if(!notPrime[i]) sushu[cnt++] = i;
for(int j=0;j<cnt&&1ll*i*sushu[j]<=n;++j) {
notPrime[i*sushu[j]] = true;
if(i%sushu[j]==0) break;
}
}
}
判素數
- 從2到sqrt(n)試除
bool solve(long long x) {
if(x<2) return false;
for(int i=2;i<=x/i;++i) {
if(x%i==0) {
return false;
}
}
return true;
}
-
大于3的素數都滿足x = 6N-1 或 x = 6N+1的性質。
我們在試除的時候用質數試除當然是更快的[N以内的素數個數大約為：N/ln(N)], 我們可以直接用6i-1和6i+1來搜，這樣的複雜度就變成了O(sqrt(N)/6)
bool isPrime(long long x) { //快速判素數
if(x<2) return false;
if(x==2||x==3) return true;
if(x%6!=5&&x%6!=1) return false;
for(int j=5;j<=x/j;j+=6) {
if(x%j==0) return false;
if(x%(j+2)==0) return false;
}
return true;
}
-
費馬素性檢驗(非确定性的)
由費馬小定理：如果P是素數，且a不是P的倍數， 那麼a^(P-1)%P = 1 ------(1)
随機找k個(k為參數)a1，a2,… ak, 帶入1式，隻要有一個不滿足，那麼就說明P不是素數。
如果都滿足，那麼有一定的機率可以認為是素數。但是存在邁克卡爾數。
不一定檢測得出來
long long fast_sum(long long a,long long b,long long mod) {
if(!a||!b) return 0;
long long ans = 0;
while(b>0) {
if(b&1) ans = (ans+a)%mod;
a = (a<<1)%mod;
}
return ans;
}
long long fast_pow(long long a,long long b,long long mod)
{
if(a==0) return 0;
long long ans = 1;
while(b>0) {
if(b&1) {
ans = fast_sum(ans,a,mod);
}
a = fast_sum(a,a,mod);
b>>=1;
}
return ans;
}
bool Fermat(long long x) { //費馬小定理：若P為素數，則a^(P-1) % P = 1 非确定性的！
if(x<2) return false;
if(x==2||x==3) return true;
int num = min(x-2,20ll);
for(int i=0;i<num;++i) {
long long a = 1ll*rand()*rand()%(x-3)+2;
if(fast_pow(a,x-1,x)!=1) {
return false;
}
}
return true;
}
-
米勒-拉賓(Miller_Rabbin)素性檢驗
基于費馬小定理 和 二次探測。基于機率的确定性的算法。複雜度很低O(log2(n))
typedef long long LL;
/*
if n < 1,373,653, it is enough to test a = 2 and 3.
if n < 9,080,191, it is enough to test a = 31 and 73.
if n < 4,759,123,141, it is enough to test a = 2, 7, and 61.
if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11.
*/
class Miller_Rabbin{
static LL fast_mul(LL a,LL b,LL mod) {
if(!a||!b) return 0;
if(mod<=2e9) { //mod<sqrt(LONG_LONG_MAX) 相乘不會溢出
if(a>=mod) a%=mod;
if(b>=mod) b%=mod;
LL res = a*b;
if(res>=mod) res%=mod;
return res;
}
LL ans = 0;
while(b>0) {
if(b&1) ans = (ans+a)%mod;
a = (a<<1)%mod;
b>>=1;
}
return ans;
}
static LL fast_pow(LL a,LL b,LL mod) {
if(a==0) return 0;
LL ans = 1;
while(b>0) {
if(b&1) ans = fast_mul(ans,a,mod);
a = fast_mul(a,a,mod);
b>>=1;
}
return ans;
}
static bool MillerRabbin(int n,int a) { //米勒拉賓素性測試
int r = 0, s = n-1, j;
if(n%a==0) return false; //a是素數，要滿足gcd(a,n)==1 隻需要滿足n%a!=0
//n和a互質
//若n是大于2的奇數，那麼n-1一定是偶數
//n-1 = 2^r * d
//費馬小定理：若n為素數，并且gcd(a,n)==1, 那麼a^(n-1) % n = 1, a^(2^r * d) % n = 1
//a^(2^r * d) = (a^d)^(2^r) = k^(2^r)
while((s&1)==0) {
s>>=1;
++r;
}
LL k = fast_pow(a,s,n); //k = a^d % n
if(k==1) return true;
for(j=0;j<r;++j,k=k*k%n) {
if(k==n-1) return true;
}
return false;
}
public:
static bool judge(int n) {
if(n<2) return false;
int su[] = {2,7,61};
for(int i=0;i<3;++i) {
if(n==su[i]) return true;
if(!MillerRabbin(n,su[i]))
return false;
}
return true;
}
};
//判斷一個數是否為素數
//測試從1到1E7
//sqrt(n)的樸素方法： 25s
//用6N+1/6N-1規律後改進的O(sqrt(n)/6) 6s
//米勒拉賓素性檢驗：3s
#include <bits/stdc++.h>
using namespace std;
const int N = 1e7+2;
typedef long long LL;
bool notPrime[N];
int sushu[700000],cnt;
//判斷一個數是否為素數
//測試從1到1E7
//sqrt(n)的樸素方法： 25s
//用6N+1/6N-1規律後改進的O(sqrt(n)/6) 6s
//米勒拉賓素性檢驗：3s
void getPrime() {
cnt = 0;
int n = N-1;
for(int i=2;i<=n;++i) {
if(!notPrime[i]) sushu[cnt++] = i;
for(int j=0;j<cnt&&1ll*i*sushu[j]<=n;++j) {
notPrime[i*sushu[j]] = true;
if(i%sushu[j]==0) break;
}
}
}
bool isPrime(LL x) { //快速判素數
if(x<2) return false;
if(x==2||x==3) return true;
if(x%6!=5&&x%6!=1) return false;
for(int j=5;j<=x/j;j+=6) {
if(x%j==0) return false;
if(x%(j+2)==0) return false;
}
return true;
}
LL fast_sum(LL a,LL b,LL mod) {
if(!a||!b) return 0;
LL ans = 0;
while(b>0) {
if(b&1) ans = (ans+a)%mod;
a = (a<<1)%mod;
}
return ans;
}
LL fast_pow(LL a,LL b,LL mod)
{
if(a==0) return 0;
LL ans = 1;
while(b>0) {
if(b&1) {
ans = fast_sum(ans,a,mod);
}
a = fast_sum(a,a,mod);
b>>=1;
}
return ans;
}
bool Fermat(LL x) { //費馬小定理：若P為素數，則a^(P-1) % P = 1
if(x<2) return false;
if(x==2||x==3) return true;
int num = min(x-2,20ll);
for(int i=0;i<num;++i) {
LL a = 1ll*rand()*rand()%(x-3)+2;
if(fast_pow(a,x-1,x)!=1) {
return false;
}
}
return true;
}
bool solve(LL x) {
if(x<2) return false;
for(int i=2;i<=x/i;++i) {
if(x%i==0) {
return false;
}
}
return true;
}
/*
if n < 1,373,653, it is enough to test a = 2 and 3.
if n < 9,080,191, it is enough to test a = 31 and 73.
if n < 4,759,123,141, it is enough to test a = 2, 7, and 61.
if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11.
*/
class Miller_Rabbin{
static LL fast_mul(LL a,LL b,LL mod) {
if(!a||!b) return 0;
if(mod<=2e9) { //mod<sqrt(LONG_LONG_MAX) 相乘不會溢出
if(a>=mod) a%=mod;
if(b>=mod) b%=mod;
LL res = a*b;
if(res>=mod) res%=mod;
return res;
}
LL ans = 0;
while(b>0) {
if(b&1) ans = (ans+a)%mod;
a = (a<<1)%mod;
b>>=1;
}
return ans;
}
static LL fast_pow(LL a,LL b,LL mod) {
if(a==0) return 0;
LL ans = 1;
while(b>0) {
if(b&1) ans = fast_mul(ans,a,mod);
a = fast_mul(a,a,mod);
b>>=1;
}
return ans;
}
static bool MillerRabbin(int n,int a) { //米勒拉賓素性測試
int r = 0, s = n-1, j;
if(n%a==0) return false; //a是素數，要滿足gcd(a,n)==1 隻需要滿足n%a!=0
//n和a互質
//若n是大于2的奇數，那麼n-1一定是偶數
//n-1 = 2^r * d
//費馬小定理：若n為素數，并且gcd(a,n)==1, 那麼a^(n-1) % n = 1, a^(2^r * d) % n = 1
//a^(2^r * d) = (a^d)^(2^r) = k^(2^r)
while((s&1)==0) {
s>>=1;
++r;
}
LL k = fast_pow(a,s,n); //k = a^d % n
if(k==1) return true;
for(j=0;j<r;++j,k=k*k%n) {
if(k==n-1) return true;
}
return false;
}
public:
static bool judge(int n) {
if(n<2) return false;
int su[] = {2,7,61};
for(int i=0;i<3;++i) {
if(n==su[i]) return true;
if(!MillerRabbin(n,su[i]))
return false;
}
return true;
}
};
int main(void) {
getPrime();
cout<<"getPrime() finish\n";
for(int i=2;i<N;++i) {
if(notPrime[i]==Miller_Rabbin::judge(i)) {
cout<<"Error "<<i<<endl;
cout<<(notPrime[i]?"no":"yes")<<endl;
return 0;
}
}
cout<<"Correct\n";
return 0;
}
關于歐拉函數
定義：歐拉函數φ(n)表示小于等于n的正整數中和n互質的數的個數。(n為正整數)
φ(1) = 1
int cal(int n) {
int ans = 0;
for(int i=1; i<=n; ++i)
if(__gcd(i,n)==1) ++ans;
return ans;
}
歐拉函數是積性函數：若m,n互質，且均為正整數，則 φ(n*m) = φ(n) * φ(m)
φ(n*m) = φ(n) * φ(m)
p為質數，則 φ(p) = p-1
φ(p) = p-1
因為如果p為質數，那麼1~p-1都和p互質，是以
φ(p) = p-1
如果n = pk, 那麼小于pk的數隻有p的倍數不和n互質。1~pk範圍内p的倍數有floor(pk/p) = pk-1個，是以φ(pk) = pk - pk-1
唯一分解定理：一個正整數n可以被唯一分解成 n = p1^a1 * p2^a2 * ... * pk^ak
n = p1^a1 * p2^a2 * ... * pk^ak
歐拉函數的求法：φ(n) = n*(1-1/p1)(1-1/p2)(1-1/p3)*(1-1/p4)……(1-1/pk)
推導過程：我們按照唯一分解定理分解n,
φ(n) = φ(p1^a1 * p2^a2 * ... * pk^ak)
因為pi都是質數，是以piai兩兩互質，由因為φ(x)是積性函數，是以
φ(n) = φ(p1^a1 * p2^a2 * ... * pk^ak) = φ(p1^a1) * φ(p2^a2) * φ(p3^a3) * ... * φ(pk^ak)
我們知道φ(pk) = pk * (1-1/p), 是以
φ(n) = φ(p1^a1 * p2^a2 * ... * pk^ak) = φ(p1^a1) * φ(p2^a2) * φ(p3^a3) * ... * φ(pk^ak) = p1^k1 * (1-1/p1) * p2^a2 * (1-1/p2) * ... * pk^ak * (1-1/pk)
=
p1^a1 * p2^a2 * ... * pk^ak * (1-1/p1) * (1-1/p2) * ... * (1-1/pk)
=
n * (1-1/p1) * (1-1/p2) * (1-1/p3) * ... * (1-1/pk)
求歐拉函數的值φ(x):
typdef long long LL;
LL Phi0(LL x) {
LL ans = x;
for(int i=2;i<=x/i;++i) {
if(x%i==0) ans = ans/i*(i-1);
while(x%i==0) x/=i;
}
if(x>1) ans = ans/x*(x-1);
return ans;
}
用6N±1優化後的代碼：
typedef long long LL;
LL Phi(LL x) {
LL ans = x;
if(x<=0) return 0;
if(x==1) return 1;
if(x%2==0) { //特判2
ans = ans/2;
while(!(x&1)) x>>=1;
}
if(x%3==0) { //特判3
ans = ans/3*2;
while(x%3==0) x/=3;
}
for(int j=5;j<=x/j;j+=6) { //用6n-1和6n+1篩
if(x%j==0) {
ans = ans/j*(j-1);
do {
x /= j;
}while(x%j==0);
}
if(x%(j+2)==0) {
ans = ans/(j+2)*(j+1);
do {
x /= j+2;
}while(x%(j+2)==0);
}
}
if(x>1) ans = ans/x*(x-1);
return ans;
}
用線性篩求歐拉函數
- 若x為質數，則Phi(x) = x-1
- 若p為質數，且i不是p的倍數，說明i和p互質，那麼
Phi(p*i) = Phi(p) * Phi(i) = (p-1)*phi[p]
- 若p為質數，且i是p的倍數，令i = pk * b (p,b互質)，那麼
Phi(i*p) = Phi(p^k*b*p) = Phi(p^(k+1)*b) = Phi(p^(k+1)) * Phi(b) = (p^(k+1)-p^k) * Phi(b) = p * (p^k-p^(k-1)) * phi(b) = p * Phi(p^k) * phi(b) = p * phi(i)
我們可以在用歐拉篩法線性篩素數的同時順便維護出歐拉函數(還可以順便維護莫比烏斯函數)
const int N = 1e7+2;
int sushu[700000],phi[N],cnt;
bool notPrime[N];
void getPhi() {
int n = N-1;
cnt = 0;
phi[1] = 1;
for(int i=2;i<=n;++i) {
if(!notPrime[i]) sushu[cnt++] = i, phi[i] = i-1;
for(int j=0;j<cnt&&1ll*i*sushu[j]<=n;++j) {
notPrime[i*sushu[j]] = true;
if(i%sushu[j]==0) {
phi[i*sushu[j]] = sushu[j]*phi[i]; //p=sushu[j] a=p^k*b,gcd(p,b)==1
//Phi(i*p) = Phi(p^(k+1)*b) = Phi(p^(k+1)*Phi(b) = (p^(k+1)-p^k)*Phi(b)
// = p*(p^k-p^(k-1))*Phi(b) = p*Phi(p^k)*Phi(b) = p*Phi(i) = phi[i]*sushu[j]
break;
} else {
phi[i*sushu[j]] = phi[i]*(sushu[j]-1); //i和sushu[j]互質，phi(sushu[j])=sushu[j]-1
}
}
}
}
求區間素數個數
鋪墊了那麼多，下面我們來求N以内素數的個數，要求N<=1E11
題目連結： LOJ-6235 區間素數個數 AC代碼<–
我們可以用min_25篩
- min_25篩是一種亞線性篩法，可以在小于O(n)的時間複雜度求出積性函數的字首和。
- 要求f( p)可以快速得到(p是素數)，f(pk)可以由f( p)快速得到
- 我們把[1,n]分成三部分：素數，合數，1
- 令prime[x]為第x個素數，prime[0]=1,prime[1]=2,prime[2]=3,prime[3]=5,prime[4]=7…,令minp(x)表示x的最小質因數
- 比如求f(x)的前n個函數值字首和，我們令g(i,j)為[2,i]範圍内f(x)的和，要求x是素數或者x的最小質因數minp(x) <= prime[j]
- 從g(i,j-1)推出g(i,j) 可以看成是埃式篩法中第j次循環，用prime[j]篩去>=prime[j]2 的合數的過程。篩去的合數有一個公因數prime[j], 我們可以提出來，1到
floor(n/prime[j])
g(i,j) = g(i,j-1)-g(n/i,i-1)
-
g(i,j) = g(i,j-1) - g(n/i,j-1) + (j-1)
-
g(i,j) = g(i,j-1)
-
g(i,j) = i-1
-
是以，我們可以二維數組遞推，也可以記憶化搜尋。
隻需要維護sqrt(n)之内的g(x,j) 和g(n/x,j)即可。
- 我們可以開兩個數組f1[], f2[], 分别記錄g(x)和g(n/x) , 也可以采用整除分塊,隻開一個g[N], 用a[]記錄整除分塊的情況。
- 第二維可以用滾動數組的方式去掉。隻需要從大到小更新g[N]即可，原理和一維的01背包相同。
我們可以把1,n的數字按照n/x的商 從小到達進行分塊，分塊的結果如下：
空間複雜度O(n0.5) 時間複雜度：O(n0.75/logn)
可以用整除分塊
#include <bits/stdc++.h>
using namespace std;
//令F(x) = (x是素數)?1:0
//那麼N以内素數的個數可以看成是f(x)的字首和
//min_25篩法可以在亞線性時間複雜度内求積性函數的字首和O(n^0.75/logn)
//在埃式篩法中，每個質數Pi從Pi*Pi,Pi*(Pi+1)一直篩到最後,每個合數必定被最小的質因數篩過一次
//我們令prime[i]為從小到大第i個素數，不妨令prime[0]=1,prime[1]=2,prime[2]=3,prime[3]=5,prime[4]=7...
//我們令g(i,j)表示前i個數,被前j個素數篩過之後剩下的數字的個數
//令minp(x)代表x的最小的質因數
//那麼g(i,j) = [1,i]中素數的個數 + [1,i]中minp()大于prime[j]的合數的個數
//g(n,0) = n-1 因為1不是素數也不是合數，先減去它，現在[2,n]都沒有别篩掉，是以g(n,0) = n-1
// { g(i,j-1) 當prime[j]*prime[j] > i的時候，因為此時prime[j]曬不掉任何i以内的數
//g(i,j) = { g(i,j-1)-g(n/prime[j],j-1)+(j-1) 埃式篩法中pj，會篩掉pj*pj,pj*(pj+1),pj*(pj+2)共有n/pj個,前面
// j-1個素數也會被重複減掉，是以再加上j-1
// { i-1 當j==0時
/*
令f(x) = 1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
g(30,0) = 29
g(2,0) = 1
g(15,0) = 14
prime[1] = 2
g(30,1) = f(2)+f(3)+f(5)+f(7)+f(9)+f(11)+...+f(29) = g(30,0)-g(15,0)+(1-1) = 29-14=15
g(10,1) = f(2)+f(3)+f(5)+f(7)+f(9) = 5 = g(10,0)-g(5,0)+(1-1) = 5
g(6,1) = f(2)+f(3)+f(5) = g(6,0)-g(6/2,0)+(1-1) = 3
g(5,1) = f(2)+f(3)+f(5) = g(5,0)-g(5/2,0)+(1-1) = 4-1 = 3
g(2,1) = f(2) = g(2,0)-g(1,0)+(1-1) = 1
prime[2] = 3
g(30,2) = f(2)+f(3)+f(5)+f(7)+f(11)+f(13)+f(17)+f(19)+f(23)+f(25)+f(29) = g(30,1)-g(30/3,1)+(2-1)
= 15-5+1 = 11 [2,3,5,7,11,13,17,19,23,25,29(11個數)]
g(6,2) = f(2)+f(3)+f(5) = g(6,1)-g(6/3,1)+(2-1) = 3
prime[3] = 5
g(30,3) = f(2)+f(3)+f(7)+f(11)+f(13)+f(17)+f(19)+f(23)+f(29) = g(30,2)-g(6,2)+(3-1)
= 11-3+2 = 10
是以30以内的質數個數為10
因為prime[3] = 5 是<= sqrt(n)的最大質數，是以後面的合數都被篩完了，可以得到答案了。
我們觀察，從g(30,1)到g(30,2)其實是埃式篩用素數prime[2]=3去篩掉了9,15,21,27這四個minp為3的數
前30個數裡面，3的倍數有[3,6,9,12,15,18,21,24,27,30] ,提取公因數後變成了
3 * [1,2,3,4,5,6,7,8,9,10], 一共有floor(30/3)個。g(10,1)=[2,3,5,7,9]=5(個),我們篩掉的其實是
3 * [2,3,5,7,9], 是[6,9,15,21,27],我們發現，在用第二個素數3去篩的時候，前面的素數會被減掉，
是以要加上(2-1)
*/
class CalPrimeCnt{
public:
static long long solve(long long n) {
int sqr = sqrt(n);
const int N = 2*sqr+10;
vector<long long> a(N);
vector<long long> g(N);
vector<int> prime(N);
//整除分塊
int m = 0; //記錄分塊的個數
for(register long long i=1;i<=n;i=a[m]+1) {
a[++m] = n/(n/i);
// printf("a[%d] = %lld\n",m,a[m]);
g[m] = a[m]-1; //g(m,0) = m-1
//a[m]記錄的是n/x得到第m大的商的最大的x值
} //分塊完了
int cnt = 0;
for(register int p=2;p<=sqr;++p) {
if(g[p]==g[p-1]) continue; //說明p為合數，f(p)=0
//用素數p去篩别人(埃式篩)
++cnt; //p為素數
long long Q = 1ll*p*p;//滾動數組，從後往前更新
//a[j]:1,2,3,4,5,6,7,10,15,30
for(register long long j=m;a[j]>=Q;--j) { //p^2=Q >= j的 g(j)不需要更新
long long block = a[j]/p; //a[j]/p = 30/2=15 n/block=30/15=2,m-2+1=m-1
//n/block = n/(a[j]/p)
int idx = block<=sqr?block:m-n/block+1;
g[j] -= g[idx]-(cnt-1); //
}
}
return g[m];
}
};
int main(void) {
// freopen("out.txt","w",stdout);
long long n = 1e11;
// cin>>n;
ios::sync_with_stdio(false);
cin.tie(0);
cout<<CalPrimeCnt::solve(n)<<"\n";
return 0;
}
（求N以内素數個數還有複雜度更優的 Meissel-Lehmer 算法，可以達到O(n(2/3)/logn) 不過據說常數略大，min_25篩也有用樹狀數組優化的O(n2/3/logn)的版本）
![POJ 1284 Primitive Roots (歐拉函數&原根定理)[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)