【BZOJ3512】DZY Loves Math IV（杜教筛）
题面
BZOJ
求
\[\sum_{i=1}^n\sum_{j=1}^m\varphi(ij)
\]
\]
其中\(n\le 10^5,m\le 10^9\)。
题解
这个数据范围很有意思。
\(n\)的值足够小，所以我们可以直接暴力枚举\(n\)。
那么所求：
\[S(n,m)=\sum_{i=1}^m\varphi(ni)
\]
\]
考虑如何将\(\varphi\)给拆开，因为\(\varphi\)只有每个质因子第一次出现的时候才会特殊计算，其他时候都直接乘起来，因此，假设\(\displaystyle n=\prod p_i^{a_i}\)，那么令\(\displaystyle p=\prod p_i^{a_i-1}\)，\(\displaystyle q=\prod p_i\)，\(pq=n\)
所以有：
\[\begin{aligned}
S(n,m)&=p\sum_{i=1}^m\varphi(qi)\\
&=p\sum_{i=1}^m \varphi(q)\varphi(\frac{i}{gcd(i,q)})gcd(i,q)\\
&=p\sum_{i=1}^m \varphi(\frac{q}{gcd(i,q)})\varphi(i)gcd(i,q)\\
&=p\sum_{i=1}^m \varphi(\frac{q}{gcd(i,q)})\varphi(i)\sum_{d|gcd(i,q)}\varphi(d)\\
&=p\sum_{i=1}^m \varphi(i)\sum_{d|i,d|q}\varphi(\frac{q}{d})\\
&=p\sum_{d|q}\varphi(\frac{q}{d})\sum_{i=1}^{m/d} \varphi(id)\\
&=p\sum_{d|q}\varphi(\frac{q}{d})S(d,[\frac{m}{d}])\\
\end{aligned}\]
S(n,m)&=p\sum_{i=1}^m\varphi(qi)\\
&=p\sum_{i=1}^m \varphi(q)\varphi(\frac{i}{gcd(i,q)})gcd(i,q)\\
&=p\sum_{i=1}^m \varphi(\frac{q}{gcd(i,q)})\varphi(i)gcd(i,q)\\
&=p\sum_{i=1}^m \varphi(\frac{q}{gcd(i,q)})\varphi(i)\sum_{d|gcd(i,q)}\varphi(d)\\
&=p\sum_{i=1}^m \varphi(i)\sum_{d|i,d|q}\varphi(\frac{q}{d})\\
&=p\sum_{d|q}\varphi(\frac{q}{d})\sum_{i=1}^{m/d} \varphi(id)\\
&=p\sum_{d|q}\varphi(\frac{q}{d})S(d,[\frac{m}{d}])\\
\end{aligned}\]
那么就可以递归处理啦？
至于复杂度？不会证啦QwQ。
#include<iostream>
#include<cstdio>
#include<vector>
#include<map>
using namespace std;
#define MOD 1000000007
#define MAX 2000000
int n,m,ans;
int pri[MAX],tot,phi[MAX],mn[MAX];
bool zs[MAX];
void pre(int n)
{
phi[1]=1;
for(int i=2;i<=n;++i)
{
if(!zs[i])pri[++tot]=i,phi[i]=i-1,mn[i]=i;
for(int j=1;i*pri[j]<=n&&j<=tot;++j)
{
zs[i*pri[j]]=true;mn[i*pri[j]]=pri[j];
if(i%pri[j])phi[i*pri[j]]=phi[i]*phi[pri[j]];
else{phi[i*pri[j]]=phi[i]*pri[j];break;}
}
}
for(int i=1;i<=n;++i)phi[i]=(phi[i-1]+phi[i])%MOD;
}
map<int,int> M,S[100010];
int Phi(int n)
{
if(n<MAX)return phi[n];
if(M[n])return M[n];
int ret=1ll*n*(n+1)/2%MOD;
for(int i=2,j;i<=n;i=j+1)
{
j=n/(n/i);
ret=(ret+MOD-1ll*(j-i+1)*Phi(n/i)%MOD)%MOD;
}
return M[n]=ret;
}
int Solve(int n,int m)
{
if(!m)return 0;
if(n==1)return Phi(m);
if(m==1)return (Phi(n)-Phi(n-1)+MOD)%MOD;
if(S[n][m])return S[n][m];
int ret=0;
vector<int> fac;int p=1,q=1,N=n;
while(N>1)
{
int x=mn[N];q*=x;N/=x;fac.push_back(x);
while(N%x==0)p*=x,N/=x;
}
for(int i=0,l=fac.size();i<1<<l;++i)
{
int d=1;
for(int j=0;j<l;++j)if(i&(1<<j))d*=fac[j];
ret=(ret+1ll*(Phi(q/d)-Phi(q/d-1)+MOD)*Solve(d,m/d))%MOD;
}
return S[n][m]=1ll*ret*p%MOD;
}
int main()
{
scanf("%d%d",&n,&m);pre(MAX);
for(int i=1;i<=n;++i)ans=(ans+Solve(i,m))%MOD;
printf("%d\n",ans);
return 0;
}