筛素数，判素数，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) 证明

p为质数，则φ(p) = p-1, φ(p^k) = p^k - p^k-1 = (p-1)*(p^k-1) = p^k * (1-1/p)

因为如果p为质数，那么1~p-1都和p互质，所以φ(p) = p-1 如果n = p^k, 那么小于p^k的数只有p的倍数不和n互质。1~p^k范围内p的倍数有floor(p^k/p) = p^k-1个，所以φ(p^k) = p^k - p^k-1

唯一分解定理：一个正整数n可以被唯一分解成n = p1^a1 * p2^a2 * ... * pk^ak的形式，其中p1,p2…pk是n的所有质因数。

欧拉函数的求法：φ(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都是质数，所以pi^ai两两互质，由因为φ(x)是积性函数，所以φ(n) = φ(p1^a1 * p2^a2 * ... * pk^ak) = φ(p1^a1) * φ(p2^a2) * φ(p3^a3) * ... * φ(pk^ak) 我们知道φ(p^k) = p^k * (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 = p^k * 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(p^k)可以由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]² 的合数的过程。筛去的合数有一个公因数prime[j], 我们可以提出来，1到floor(n/prime[j])即为筛去的数的范围。g(i,j) = g(i,j-1)-g(n/i,i-1),但是这样的话，g(n/i,j-1)中也包含了前j-1个素数, 我们多减了，所以药加回来。还要加上(j-1)

g(i,j) = g(i,j-1) - g(n/i,j-1) + (j-1) （当prime[j]² <= n的时候）

g(i,j) = g(i,j-1) （当prime[j]² > n的时候，埃式筛法的第j次循环没有筛掉任何一个合数）

g(i,j) = i-1 （当j==0的时候）

所以，我们可以二维数组递推，也可以记忆化搜索。 只需要维护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(n^0.5) 时间复杂度：O(n^0.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(n^2/3/logn)的版本）

---