[BZOJ 3512] DZY Loves Math IV
一、题目
给定 \(n,m\) ，求下面的柿子模 \(\tt 1e9+7\) 的值：
\[\sum_{i=1}^n\sum_{j=1}^m\varphi(ij)
\]
\(1\leq n\leq1e5,1\leq m\leq 1e9\)
二、解法
发现 \(n\) 很小，可以尝试枚举 \(n\) 这一维，加速 \(m\) 那一位的运算，令 \(S(n,m)=\sum_{i=1}^m\varphi(ni)\) ：
\[\sum_{i=1}^m\varphi(ni)
\]
如果我们想用杜教筛快速求的话，是必须要把这个东西拆开的，利用 \(\varphi\) 是积性函数，可以知道我们必须要把他们转化成互质的两个数才能拆，设 \(w=\prod p_i\)，也就是 \(n\) 质因数分解的质数一次乘积，\(y=\frac{n}{w}\) ，则：
\[=\sum_{i=1}^n\varphi(wyi)=y\sum_{i=1}^n\varphi(wi)
\]
但是现在好像还是不能拆，我们尝试引入 \(\gcd(w,i)\) ，这样也许就可以互质了：
\[=y\sum_{i=1}^n\varphi(\frac{w}{\gcd(w,i)}\times(i\times gcd(w,i)))=y\sum_{i=1}^n\varphi(\frac{w}{\gcd(w,i)})\varphi(i\times\gcd(w,i))
\]
\[=y\sum_{i=1}^n\varphi(\frac{w}{\gcd(w,i)})\varphi(i)\gcd(w,i)
\]
看到 \(\gcd\) 就有点想反演他，欧拉反演出奇迹：
\[=y\sum_{i=1}^n\varphi(\frac{w}{\gcd(w,i)})\varphi(i)\sum_{e|\gcd(w,i)}\varphi(e)
\]
因为 \(\frac{w}{\gcd(w,i)}\) 和 \(\gcd(w,i)\) 互质，所以他和 \(e\) 也互质，我们把他们结合：
\[=y\sum_{i=1}^n\varphi(i)\sum_{e|\gcd(w,i)}\varphi(\frac{w}{\frac{\gcd(w,i)}{e}})=y\sum_{i=1}^n\varphi(i)\sum_{e|\gcd(w,i)}\varphi(\frac{w}{e})
\]
可以先枚举 \(e\) 了，根据莫比乌斯反演的经验这样做一定有用：
\[=y\sum_{e|n}\varphi(\frac{w}{e})\sum_{i=1}^{m/e}\varphi(ie)=y\sum_{e|n}\varphi(\frac{w}{e})S(e,m/e)
\]
这样就推出了 一个递归的问题 ，递归的边界是 \(n=1\) ，可以用杜教筛解决的。而且由于第二维的所有取值一定是某个 \(\frac{m}{x}\) ，所以杜教筛只会完整地筛一次，时间复杂度 \(O(m^{\frac{2}{3}})\) 。然后对于算 \(S\) 的复杂度由于第一维一定是质数的一次乘积（一开始传进去的 \(n\) 单独判一下就行了），所以也不会很多，第二维的取值只有 \(\sqrt m\) 个。
#include <cstdio>
#include <iostream>
#include <vector>
#include <map>
using namespace std;
const int N = 1000000;
const int M = N+5;
const int MOD = 1e9+7;
#define int long long
int read()
{
int x=0,f=1;char c;
while((c=getchar())<'0' || c>'9') {if(c=='-') f=-1;}
while(c>='0' && c<='9') {x=(x<<3)+(x<<1)+(c^48);c=getchar();}
return x*f;
}
int n,m,cnt,ans,p[M],sp[M],phi[M],low[M],nmsl;
map<int,int> sum,zxy[M];
void init(int n)
{
phi[1]=low[1]=1;
for(int i=2;i<=n;i++)
{
if(!phi[i])
{
low[i]=i;//表示i的一次质数乘积
phi[i]=i-1;
p[++cnt]=i;
}
for(int j=1;j<=cnt && i*p[j]<=n;j++)
{
low[i*p[j]]=i;
if(i%p[j]==0)
{
phi[i*p[j]]=phi[i]*p[j];
low[i*p[j]]=low[i];
break;
}
phi[i*p[j]]=phi[i]*(p[j]-1);
low[i*p[j]]=low[i]*p[j];
}
}
for(int i=1;i<=n;i++)
sp[i]=(phi[i]+sp[i-1])%MOD;
}
int zy(int x)
{
return x*(x+1)/2%MOD;
}
int get(int n)
{
if(n<=N) return sp[n];
if(sum[n]) return sum[n];
int ans=zy(n);
for(int l=2,r;l<=n;l=r+1)
{
r=n/(n/l);
ans=(ans-(r-l+1)*get(n/l))%MOD;
}
return sum[n]=ans;
}
int jzm(int n,int m)
{
if(n==1) return get(m);
if(m==1) return phi[n];
if(!m) return 0;
if(zxy[n][m]) return zxy[n][m];
int r=0;
for(int i=1;i*i<=n;i++)
if(n%i==0)
{
r=(r+phi[n/i]*jzm(i,m/i))%MOD;
if(i*i!=n) r=(r+phi[i]*jzm(n/i,m/(n/i)));
}
return zxy[n][m]=r;
}
signed main()
{
init(N);
n=read();m=read();
for(int i=1;i<=n;i++)
ans=(ans+(i/low[i])*jzm(low[i],m))%MOD;
printf("%lld\n",(ans+MOD)%MOD);
}
