杜教筛
运用狄利克雷卷积的形式对一些积性函数在小于线性的时间内求前缀和。
套路式
对于要求前缀和的积性函数(f(i))，假设其前缀和函数为(S(i))。构造积性函数(g(i))，与原函数做狄利克雷卷积得$$(f*g)(i)=sum_{d|i}g(d)f(frac id)$$
对卷积函数求前缀和
[sum_{i=1}^{n}(f*g)(i)=sum_{i=1}^{n}sum_{d|i}g(d)f(frac id)\=sum_{d=1}^{n}g(d)sum_{d|i}^{n}f(frac id)\=sum_{d=1}^{n}g(d)sum_{i=1}^{n/d}f(i)\=sum_{d=1}^{n}g(d)S(frac nd)
]
所以就得到这样一个等式
[g(1)S(n)=sum_{i=1}^{n}(f*g)(i)-sum_{d=2}^{n}g(d)S(frac nd)
]
e.g.
求(sum_{i=1}^{n}varphi(i),nle10^{10})
设(S(x)=sum_{i=1}^{x}varphi(i))，则有
[g(1)S(n)=sum_{i=1}^{n}(f*g)(i)-sum_{d=2}^{n}g(d)S(frac nd)
]
我们要构造一个函数使得(sum_{i=1}^{n}(f*g)(i))很好求
利用(sum_{d|n}varphi(d)=n)可知令(g(i)=1)就可以有
[(f*g)(i)=sum_{d|i}varphi(frac id)=sum_{d|i}varphi(d)=i
]
[sum_{i=1}^{n}(f*g)(i)=sum_{i=1}^{n}i=frac{i(i+1)}{2}
]
所以就搞定了。
关于实现
写一个求(S(i))的函数，计算的时候递归调用。注意记忆化。
一般而言都是开map
Gay神说可以不开map省掉一个log
以下出自Gay神：
(qquad)杜教筛其实有个小trick,可以不用写map,省掉一个log,具体是因为你所有要筛出来的前缀和都是可以写成F(n/i)的,前面预处理了的都不用记了,你对于i记忆化,就可以达到快速记下要筛出来的前缀和的目的.
所以开一个数组记录就可以了。
（注：经yyb实测表明这并不能优化多少，所以以下两份代码均未使用这种trick）
模板题
【BZOJ3944】Sum
就是求上面讲到的两个东西
#include<cstdio>
#include<algorithm>
#include<map>
using namespace std;
#define ll long long
const int N = 2500000;
ll gi()
{
ll x=0,w=1;char ch=getchar();
while ((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
if (ch=='-') w=0,ch=getchar();
while (ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
return w?x:-x;
}
int pri[N+5],tot,zhi[N+5];
ll mu[N+5],phi[N+5];
map<ll,ll>M,P;
void Mobius()
{
zhi[1]=1;mu[1]=phi[1]=1;
for (int i=2;i<=N;++i)
{
if (!zhi[i]) pri[++tot]=i,mu[i]=-1,phi[i]=i-1;
for (int j=1;j<=tot&&i*pri[j]<=N;++j)
{
zhi[i*pri[j]]=1;
if (i%pri[j]) mu[i*pri[j]]=-mu[i],phi[i*pri[j]]=phi[i]*phi[pri[j]];
else {mu[i*pri[j]]=0;phi[i*pri[j]]=phi[i]*pri[j];break;}
}
}
for (int i=1;i<=N;++i)
mu[i]+=mu[i-1],phi[i]+=phi[i-1];
}
ll Phi(ll x)
{
if (x<=N) return phi[x];
if (P[x]) return P[x];
ll res=0,i=2,j;
while (i<=x)
{
j=x/(x/i);
res+=(j-i+1)*Phi(x/i);
i=j+1;
}
return P[x]=x*(x+1)/2-res;
}
ll Mu(ll x)
{
if (x<=N) return mu[x];
if (M[x]) return M[x];
ll res=0,i=2,j;
while (i<=x)
{
j=x/(x/i);
res+=(j-i+1)*Mu(x/i);
i=j+1;
}
return M[x]=1-res;
}
int main()
{
Mobius();
int T=gi();
while (T--)
{
int n=gi();
printf("%lld %lld
",Phi(n),Mu(n));
}
return 0;
}
还是模板题
【BZOJ4916】神犇和蒟蒻
求(sum_{i=1}^{n}mu(i^2))和(sum_{i=1}^{n}varphi(i^2))
第一问答案恒为1.
第二问，(varphi(i^2)=ivarphi(i))，所以还是同样的套路做。
(sum_{i=1}^{n}i^2=frac{1}{6}n(n+1)(2n+1))
#include<cstdio>
#include<algorithm>
#include<map>
using namespace std;
const int mod = 1e9 + 7;
const int N = 10000000;
const int inv = 166666668;
int gi()
{
int x=0,w=1;char ch=getchar();
while ((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
if (ch=='-') w=0,ch=getchar();
while (ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
return w?x:-x;
}
int pri[N+5],tot,zhi[N+5],phi[N+5];
map<int,int>P;
void Mobius()
{
zhi[1]=phi[1]=1;
for (int i=2;i<=N;i++)
{
if (!zhi[i]) pri[++tot]=i,phi[i]=i-1;
for (int j=1;j<=tot&&i*pri[j]<=N;j++)
{
zhi[i*pri[j]]=1;
if (i%pri[j]) phi[i*pri[j]]=1ll*phi[i]*phi[pri[j]]%mod;
else {phi[i*pri[j]]=1ll*phi[i]*pri[j]%mod;break;}
}
}
for (int i=1;i<=N;i++) phi[i]=(phi[i-1]+1ll*i*phi[i]%mod)%mod;
}
int Phi(int n)
{
if (n<=N) return phi[n];
if (P[n]) return P[n];
int res=1ll*n*(n+1)%mod*(2*n+1)%mod*inv%mod;
int i=2,j;
while (i<=n)
{
j=n/(n/i);
res=(res-1ll*(j+i)*(j-i+1)/2%mod*Phi(n/i)+mod)%mod;
i=j+1;
}
return P[n]=res;
}
int main()
{
Mobius();
int n=gi();
printf("1
%d
",Phi(n));
return 0;
}