模板 - 杜教筛
取模版本:
const int mod=1e9+7;
const int inv2=(mod+1)>>1;
const int MAXN=5e6;
unordered_map<ll,ll> Smu;
unordered_map<ll,ll> Sphi;
ll sphi[MAXN+1];
ll smu[MAXN+1];
int pri[MAXN+1];
int &pritop=pri[0];
void sieve(int n=MAXN) {
sphi[1]=1;
smu[1]=1;
for(int i=2; i<=n; i++) {
if(!pri[i]) {
pri[++pritop]=i;
sphi[i]=i-1;
smu[i]=-1;
}
for(int j=1; j<=pritop; j++) {
int &p=pri[j];
int t=ip;
if(t>n)
break;
pri[t]=1;
if(i%p) {
sphi[t]=sphi[i]sphi[p];
if(sphi[t]>=mod)
sphi[t]%=mod;
smu[t]=-smu[i];
} else {
sphi[t]=sphi[i]*p;
if(sphi[t]>=mod)
sphi[t]%=mod;
smu[t]=0;
break;
}
}
}
for(int i=2; i<=n; i++) {
sphi[i]+=sphi[i-1];
if(sphi[i]>=mod)
sphi[i]%=mod;
smu[i]+=smu[i-1];
if(smu[i]<0)
smu[i]+=mod;
}
}
inline ll s1(ll n) {
if(n>=mod)
n%=mod;
ll tmp=n(n+1);
if(tmp>=mod)
tmp%=mod;
tmp=inv2;
if(tmp>=mod)
tmp%=mod;
return tmp;
}
inline ll Mu(ll n) {
if(n<=MAXN)
return smu[n];
if(Smu.count(n))
return Smu[n];
ll ret=1;
for(ll l=2,r; l<=n; l=r+1) {
ll t=n/l;
r=n/t;
ll tmp=r-(l-1);
if(tmp<0)
tmp+=mod;
tmp*=Mu(t);
if(tmp>=mod)
tmp%=mod;
ret-=tmp;
if(ret<0)
ret+=mod;
}
return Smu[n]=ret;
}
inline ll Phi(ll n) {
if(n<=MAXN)
return sphi[n];
if(Sphi.count(n))
return Sphi[n];
ll ret=s1(n);
for(ll l=2,r; l<=n; l=r+1) {
ll t=n/l;
r=n/t;
ll tmp=r-(l-1);
tmp*=Phi(t);
if(tmp>=mod)
tmp%=mod;
ret-=tmp;
if(ret<0)
ret+=mod;
}
return Sphi[n]=ret;
}
int main() {
sieve();
ll n;
scanf("%lld",&n);
printf("%lld\n",Phi(n));
}
</details>
---
2019/6/2更新,杜教筛可以省去哈希表的复杂度和空间。当n不改变的时候,可以发现我们关心的都是n/i,大于sqrt(n)的是不超过sqrt(n)个。
可以换一种记忆方式省去哈希表的复杂度。(虽然实际上好像并没多快,反而还慢了,还巨难写?)但是当n改变时,需要清除多余的记忆。
<details>
```cpp
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int INF=0x3f3f3f3f;
//N为n^(2/3)最快
int n;
const int N23=pow(1.0*INT_MAX,2.0/3.0);
const int N12=pow(1.0*INT_MAX,1.0/2.0);
const int MAXN=N23;
ll sump[N23+N12+5];
int sumu[N23+N12+5];
int pri[MAXN+5],pritop;
bool notpri[MAXN+5];
//pritop从1开始计数
void sieve(int n=MAXN) {
//对于欧拉函数,-1是特殊值,对于莫比乌斯函数,INF是特殊值
notpri[0]=1,sump[0]=0,sumu[0]=0;
notpri[1]=1,sump[1]=1,sumu[1]=1;
for(int i=2; i<=n; i++) {
if(!notpri[i])
pri[++pritop]=i,sump[i]=i-1,sumu[i]=-1;
for(int j=1; j<=pritop&&i*pri[j]<=n; j++) {
notpri[i*pri[j]]=1;
//略有不同
if(i%pri[j])
sump[i*pri[j]]=sump[i]*sump[pri[j]],sumu[i*pri[j]]=-sumu[i];
else {
sump[i*pri[j]]=sump[i]*pri[j],sumu[i*pri[j]]=0;;
break;
}
}
}
for(int i=2; i<=n; i++) {
sump[i]+=sump[i-1];
sumu[i]+=sumu[i-1];
}
for(int i=MAXN+1;i<=MAXN+N12;i++){
sump[i]=-1;
sumu[i]=INF;
}
}
inline int idx(int x) {
return (x<=MAXN)?x:(n/x+MAXN);
}
//杜教筛欧拉函数
inline ll GetSphi(int n) {
int t=idx(n);
if(sump[t]!=-1)
return sump[t];
//小心溢出
ll ret=((1ll+n)*n)/2;
//f*g=id的前缀和
for(int l=2,r; r<2147483647&&l<=n; l=r+1) {
r=n/(n/l);
ret-=(r-l+1)*GetSphi(n/l);
//同上,因为两个的g都是I
}
return sump[t]=ret; // 记忆化
}
//杜教筛莫比乌斯函数
inline int GetSumu(int n) {
int t=idx(n);
if(sumu[t]!=INF)
return sumu[t];
int ret=1; // 单位元的前缀和就是 1
for(int l=2,r;r<2147483647&&l<=n;l=r+1) {
r=n/(n/l);
ret-=(r-l+1)*GetSumu(n/l);
//(r-l+1)就是I在[l,r]的和
}
return sumu[t]=ret;// 记忆化
}
int main() {
#ifdef Yinku
freopen("Yinku.in","r",stdin);
#endif // Yinku
sieve();
int t;
scanf("%d",&t);
while(t--) {
for(int i=MAXN+1;i<=MAXN+N12;i++){
sump[i]=-1;
sumu[i]=INF;
}
scanf("%d",&n);
printf("%lld %d\n",GetSphi(n),GetSumu(n));
}
}
因为不会杜教筛(还有质数间隔),被一万六学院调教了。
补充一波杜教筛。
#include<bits/stdc++.h>
#define ll long long
using namespace std;
//N为n^(2/3)最快
int n;
const int MAXN=5e6;
unordered_map<int,int> Smu;
unordered_map<int,ll> Sphi;
ll sump[MAXN+5];
int sumu[MAXN+5];
int pri[MAXN+5],pritop;
bool notpri[MAXN+5];
//pritop从1开始计数
void sieve(int n) {
notpri[1]=sump[1]=sumu[1]=1;
for(int i=2; i<=n; i++) {
if(!notpri[i])
pri[++pritop]=i,sump[i]=i-1,sumu[i]=-1;
for(int j=1; j<=pritop&&i*pri[j]<=n; j++) {
notpri[i*pri[j]]=1;
//略有不同
if(i%pri[j])
sump[i*pri[j]]=sump[i]*sump[pri[j]],sumu[i*pri[j]]=-sumu[i];
else {
sump[i*pri[j]]=sump[i]*pri[j],sumu[i*pri[j]]=0;;
break;
}
}
}
for(int i=2; i<=n; i++) {
sump[i]+=sump[i-1];
sumu[i]+=sumu[i-1];
}
}
//杜教筛莫比乌斯函数
inline ll GetSumu(int n) {
if(n <= MAXN)
return sumu[n]; // sumu是提前筛好的前缀和
if(Smu.count(n))
return Smu[n]; // 记忆化
ll ret = 1ll; // 单位元的前缀和就是 1
for(int l = 2, r; r<2147483647&&l <= n; l = r + 1) {
r = n / (n / l);
ret -= (r - l + 1) * GetSumu(n / l);
// (r - l + 1) 就是 I 在 [l, r] 的和
}
return Smu[n] = ret; // 记忆化
}
//杜教筛欧拉函数
inline ll GetSphi(int n) {
if(n <= MAXN)
return sump[n]; // 提前筛好的
if(Sphi.count(n)) return Sphi[n]; // 记忆化
ll ret = n;
if(ret%2==0)
ret=(ret/2)* (1ll+n);
else
ret*=((1ll+n)/2);
// f * g = id 的前缀和
for(int l = 2, r; r<2147483647&&l <= n; l = r + 1) {
r = n / (n / l);
ret -= (r - l + 1) * GetSphi(n / l);
// 同上,因为两个的 g 都是 I
}
return Sphi[n] = ret; // 记忆化
}
int main() {
sieve(MAXN);
int t;
scanf("%d",&t);
while(t--) {
scanf("%d",&n);
printf("%lld %lld\n",GetSphi(n),GetSumu(n));
}
}
用来快速求出积性函数的前缀和(因为一般分块也是用前缀和,所以很多莫比乌斯反演可以用杜教筛优化到 \(O(n^{\frac{2}{3}})\) 预处理,然后整除分块 \(O(n^{\frac{1}{2}})\) 查询)。
前置知识
积性函数
积性函数:对于任意互质的整数 \(a,b\) 有 \(f(ab)=f(a)f(b)\) 则称 \(f(x)\) 的数论函数。
完全积性函数:对于任意整数 \(a,b\) 有 \(f(ab)=f(a)f(b)\) 的数论函数。
积性函数的积未必是积性函数,狄利克雷卷积才是。
常见的积性函数:\(\varphi,\mu,\sigma,d\)
常见的完全积性函数:\(\epsilon,I,id\)
这里特殊解释一下 \(\epsilon,I,id\) 分别是什么意思: \(\epsilon(n) = [n=1], I(n) = 1, id(n) = n\)
狄利克雷卷积
设 \(f, g\) 是两个数论函数,它们的狄利克雷卷积卷积是:\((f*g)(n) = \sum \limits _{d | n} f(d) g(\frac{n}{d})\)
性质:满足交换律,结合律
单位元: \(\epsilon\) (即 \(f*\epsilon=f\) )
结合狄利克雷卷积得到的几个性质:
\(\mu * I = \epsilon\)
\(\varphi * I = id\)
\(\mu * id = \varphi\)
杜教筛
设现在要求积性函数 \(f\) 的前缀和, 设 \(\sum \limits_{i=1}^{n} f(i) = S(n)\) 。
再找一个积性函数 \(g\) ,则考虑它们的狄利克雷卷积的前缀和
\(\sum \limits_{i=1}^{n}(f*g)(i)\)
\(= \sum \limits_{i=1}^{n} \sum \limits_{d|i} f(d)g(\frac{i}{d})\)
\(= \sum \limits_{d=1}^{n} g(d)\sum\limits_{i=1}^{\lfloor \frac{n}{d}\rfloor } f(i)\)
\(= \sum \limits_{d=1}^{n} g(d) S(\lfloor \frac{n}{d} \rfloor)\)
则:
\(g(1)S(n)=\sum\limits_{i=1}^{n}(f*g)(i) - \sum \limits _{i=2}^{n} g(i) S(\lfloor \frac{n}{i} \rfloor)\)
//杜教筛莫比乌斯函数
inline ll GetSumu(int n) {
if(n <= N) return sumu[n]; // sumu是提前筛好的前缀和
if(Smu[n]) return Smu[n]; // 记忆化
ll ret = 1ll; // 单位元的前缀和就是 1
for(int l = 2, r; l <= n; l = r + 1) {
r = n / (n / l); ret -= (r - l + 1) * GetSumu(n / l);
// (r - l + 1) 就是 I 在 [l, r] 的和
} return Smu[n] = ret; // 记忆化
}
//杜教筛欧拉函数
inline ll GetSphi(int n) {
if(n <= N) return sump[n]; // 提前筛好的
if(Sphi[n]) return Sphi[n]; // 记忆化
ll ret = 1ll * n * (n + 1) / 2; // f * g = id 的前缀和
for(int l = 2, r; l <= n; l = r + 1) {
r = n / (n / l); ret -= (r - l + 1) * GetSphi(n / l);
// 同上,因为两个的 g 都是 I
} return Sphi[n] = ret; // 记忆化
}
注意一些常用的变换,例如 \(\sum\limits_{d|n}\mu(n)=[n==1]\) ,\(\sum\limits_{d|n}\varphi(n)=n\)
不用哈希表版本还带取模的杜教筛经过了欧拉函数前缀和的验证,跑得飞快。
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int INF=0x3f3f3f3f;
const int mod=1e9+7;
const int inv2=(mod+1)>>1;
//N为n^(2/3)最快
ll n;
const int MAXN23=pow(1e10,2.0/3.0);
const int MAXN12=pow(1e10,1.0/2.0);
const int MAXN=MAXN23;
ll sump[MAXN23+MAXN12+10];
int pri[MAXN+10],pritop;
bool notpri[MAXN+10];
//pritop从1开始计数
void sieve(int n=MAXN) {
notpri[0]=1,sump[0]=0;
notpri[1]=1,sump[1]=1;
for(int i=2; i<=n; i++) {
if(!notpri[i])
pri[++pritop]=i,sump[i]=(i-1);
for(int j=1,t; j<=pritop&&(t=i*pri[j])<=n; j++) {
notpri[t]=1;
//略有不同
if(i%pri[j]){
sump[t]=sump[i]*sump[pri[j]];
if(sump[t]>=mod)
sump[t]%=mod;
}
else {
sump[t]=sump[i]*pri[j];
if(sump[t]>=mod)
sump[t]%=mod;
break;
}
}
}
//cout<<"done!"<<endl;
for(int i=2; i<=n; i++) {
sump[i]+=sump[i-1];
if(sump[i]>=mod)
sump[i]-=mod;
}
//不用考虑边界,因为比MAXN还大的才会进入这个记忆化
for(int i=MAXN+1;i<=MAXN+MAXN12;i++){
sump[i]=-1;
}
}
inline int idx(ll x) {
return (x<=MAXN)?x:(n/x+MAXN);
}
inline ll s1(ll n){
if(n>=mod)
n%=mod;
ll t=(n*(n+1)/2);
if(t>=mod)
t%=mod;
return t;
}
//杜教筛欧拉函数
inline ll GetSphi(ll n) {
int t=idx(n);
if(sump[t]!=-1)
return sump[t];
//小心溢出
ll ret=s1(n);
//f*g=id的前缀和
for(ll l=2,r;l<=n; l=r+1) {
r=n/(n/l);
ll tmp=(r-l+1)%mod*GetSphi(n/l)%mod;
ret-=tmp;
if(ret<0)
ret+=mod;
//同上,因为两个的g都是I
}
return sump[t]=ret; // 记忆化
}
int main() {
#ifdef Yinku
freopen("Yinku.in","r",stdin);
#endif // Yinku
sieve();
scanf("%lld",&n);
printf("%lld",GetSphi(n));
}