Miller_Rabin 判素数
Pollard_Rho 分解质因子
区间筛
- 筛出某区间 ([L,R]) 内的所有质数/所有数的某积性函数。 (L,Rleq 10^{12},R-Lleq 10^6)
预处理区间 ([1,10^6]) 内的质数。对 ([1,10^6]) 内任意质数 (p) ，将区间 ([L,R]) 内 (p) 的倍数，按埃氏筛的方式筛去/算出该质因子对积性函数的贡献。最后遍历区间 ([L,R]) ，找出区间质数/算出剩余质因子对积性函数的贡献。
例题 1：洛谷 P1835 素数密度 求区间质数：预处理 ([1,10^6]) 内质数，筛去 ([L,R]) 内非质数。
#include<bits/stdc++.h>
using namespace std;
const int MAXN=1e6+10;
int L,R,isnotprime[MAXN],Isnotprime[MAXN];
inline void sieve(){
Isnotprime[1]=1;
for(int i=2;i<=1e6;i++)
if(!Isnotprime[i]){
for(int j=i+i;j<=1e6;j+=i)
Isnotprime[j]=1;
long long num=L+(i-L%i)%i,pos;
if(num==i) num+=i;
for(pos=num-L+1;num<=R;num+=i,pos+=i)
isnotprime[pos]=1;
}
}
int main(){
cin>>L>>R;
sieve();
int Cnt=0;
for(long long i=1,j=L;j<=R;i++,j++)
if(!isnotprime[i]) Cnt++;
cout<<Cnt;
return 0;
}
例题 2：洛谷 P3601 签到题 求区间 (displaystyle sum_{i=l}^r[i-oldsymbolvarphi(i)]=sum(r)-sum(l-1)-sum_{i=l}^roldsymbol varphi(i),sum(n)={n(n+1)over 2})
等价于求区间欧拉函数：预处理 ([1,10^6]) 内质数，筛去 ([L,R]) 内质因子小于 (10^6) 的质因子部分。最后的遍历，对剩余的质因子算出对欧拉函数的贡献
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAXN=1e6+10;
const ll MOD=666623333;
ll L,R,Phi[MAXN],Val[MAXN];
bool Isnotprime[MAXN];
inline void sieve(){
Isnotprime[1]=1;
for(int i=2;i<=1e6;i++) if(!Isnotprime[i]){
for(int j=i+i;j<=1e6;j+=i) Isnotprime[j]=1;
ll num=L+(i-L%i)%i,pos;
if(num==i) num+=i;
for(pos=num-L+1;num<=R;num+=i,pos+=i){
Phi[pos]=Phi[pos]/i*(i-1);
while(Val[pos]%i==0) Val[pos]/=i;
}
}
for(ll i=1,j=L;j<=R;i++,j++)
if(Val[i]!=1) Phi[i]=Phi[i]/Val[i]*(Val[i]-1);//剩余质因子不为 1
}
inline __int128 sum(__int128 n){
return n*(n+1)/2;
}
int main(){
cin>>L>>R;
for(ll i=1,j=L;j<=R;i++,j++) Phi[i]=Val[i]=j;
sieve();
ll Ans=(sum(R)-sum(L-1))%MOD;
for(ll i=1,j=L;j<=R;i++,j++) Ans=(Ans-Phi[i])%MOD;
Ans=(Ans+MOD)%MOD;
cout<<Ans;
return 0;
}
- 筛出某区间 ([L,R]) 内的所有质数/所有数的某积性函数。 (L,Rleq 10^{18},R-Lleq 10^6)
预处理区间 ([1,10^6]) 内的质数。对 ([1,10^6]) 内任意质数 (p) ，将区间 ([L,R]) 内 (p) 的倍数，按埃氏筛的方式筛去/算出该质因子对积性函数的贡献。剩余的质因子形式有三种：(p,p^2,pq(p,qin Primewedge p,qgeq 10^6)) 。
最后遍历区间 ([L,R]) ，对未筛去的数，先直接开方，判定是否为 (p^2) 形式；再用 Miller_Rabin 算法判质数/用 Pollard_Rho 算法求出最小质因子，进而求出对积性函数的贡献。
例题：洛谷 P3653 小清新数学题 求区间莫比乌斯函数：先筛出 ([1,10^6]) 内素数，再算出区间 ([L,R]) 内小质因子对莫比乌斯函数的贡献。最后遍历 ([L,R]) 区间，先直接开方，判定是否为 (p^2) 形式，是则莫比乌斯函数直接为 (0) ；再用 Miller_Rabin 算法判质数，是则莫比乌斯函数变为相反数，否则为 (pq) 形式，(oldsymbolmu(p)oldsymbolmu(q)=1) ，莫比乌斯函数值不变
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll MAXN=1e6+10;
ll L,R,Mu[MAXN],Val[MAXN];
bool isnotprime[MAXN];
inline ll fpow(ll a,ll x,ll mod) { ll ans=1; for(;x;x>>=1,a=(__int128)a*a%mod) if(x&1) ans=(__int128)ans*a%mod; return ans; }
//快速幂
inline ll lowbit(ll val) { return val&(-val); }
inline bool Miller_Rabin(ll val,ll p){
if( fpow(p,val-1,val)!=1 ) return 0;
ll lim=lowbit(val-1),num=fpow(p,(val-1)/lim,val),Sqr;
for(ll bit=1;bit<=lim;bit<<=1){
Sqr=(__int128)num*num%val;
if(Sqr==1&&num!=1&&num!=val-1) return 0;
num=Sqr;
}
return 1;
}
inline bool isprime(ll val){
for(int i=1;i<=37;i++) if(!isnotprime[i]){
if(val%i==0) return val==i;
if(!Miller_Rabin(val,i)) return 0;
}
return 1;
}
inline void sieve(){
isnotprime[1]=1;
for(ll num=L,pos=1;num<=R;num++,pos++) Mu[pos]=1,Val[pos]=num;
for(int i=2;i<=1e6;i++) if(!isnotprime[i]){
for(int j=i+i;j<=1e6;j+=i) isnotprime[j]=1;
ll num=L+(i-L%i)%i,pos;
if(num==i) num+=i;
for(pos=num-L+1;num<=R;pos+=i,num+=i){
int Cnt=0;
while(Val[pos]%i==0) Val[pos]/=i,Cnt++;
if(Cnt==1) Mu[pos]=-Mu[pos];
else Mu[pos]=0;
}
}
for(ll num=L,pos=1;num<=R;num++,pos++) if(Val[pos]!=1) {
ll tmp=sqrt(Val[pos])+0.5;
if(tmp*tmp==Val[pos]) Mu[pos]=0;//p^2 形式
else if( isprime(Val[pos]) ) Mu[pos]=-Mu[pos];//p 形式
//pq 形式不变
}
}
int main(){
cin>>L>>R;
sieve();
ll Ans=0;
for(ll num=L,pos=1;num<=R;num++,pos++) Ans+=Mu[pos];
cout<<Ans;
return 0;
}
区间筛的一般形式：
inline bool Miller_Rabin(ll val,ll p) {...}
inline bool isprime(ll val) {...}
inline ll Pollar_Rho(ll val) {...}
inline ll f(ll p,ll k) {...}//求解 p^k 积性函数
bool isnotprime[MAXN];
ll F[MAXN],Val[MAXN];
inline void sieve(){
isnotprime[1]=1;
for(ll num=L,pos=1;num<=R;num++,pos++) F[pos]=1,Val[pos]=num;
for(int i=2;i<=1e6;i++) if(!isnotprime[i]){
for(int j=i+i;j<=1e6;j+=i) isnotprime[j]=1;
ll num=L+(i-L%i)%i,pos;
if(num==i) num+=i;
for(pos=num-L+1;num<=R;pos+=i,num+=i){
int Cnt=0;
while(Val[pos]%i==0) Val[pos]/=i,Cnt++;
F[pos]*=f(i,Cnt);
}
}
for(ll num=L,pos=1;num<=R;num++,pos++) if(Val[pos]!=1) {
ll tmp=sqrt(Val[pos])+0.5;
if(tmp*tmp==Val[pos]) F[pos]*=f(tmp,2);//p^2 形式
else if( isprime(Val[pos]) ) F[pos]*=f(Val[pos],1);//p 形式
else{//pq 形式不变
tmp=Pollard_Rho(Val[pos]);
F[pos]*=f(tmp,1);
F[pos]*=f(Val[pos]/tmp,1);
}
}
}