一、题目
给定 (n,m) ，求下面的柿子模 ( t 1e9+7) 的值：
[sum_{i=1}^nsum_{j=1}^mvarphi(ij)
]
(1leq nleq1e5,1leq mleq 1e9)
二、解法
发现 (n) 很小，可以尝试枚举 (n) 这一维，加速 (m) 那一位的运算，令 (S(n,m)=sum_{i=1}^mvarphi(ni)) ：
[sum_{i=1}^mvarphi(ni)
]
如果我们想用杜教筛快速求的话，是必须要把这个东西拆开的，利用 (varphi) 是积性函数，可以知道我们必须要把他们转化成互质的两个数才能拆，设 (w=prod p_i)，也就是 (n) 质因数分解的质数一次乘积，(y=frac{n}{w}) ，则：
[=sum_{i=1}^nvarphi(wyi)=ysum_{i=1}^nvarphi(wi)
]
但是现在好像还是不能拆，我们尝试引入 (gcd(w,i)) ，这样也许就可以互质了：
[=ysum_{i=1}^nvarphi(frac{w}{gcd(w,i)} imes(i imes gcd(w,i)))=ysum_{i=1}^nvarphi(frac{w}{gcd(w,i)})varphi(i imesgcd(w,i))
]
[=ysum_{i=1}^nvarphi(frac{w}{gcd(w,i)})varphi(i)gcd(w,i)
]
看到 (gcd) 就有点想反演他，欧拉反演出奇迹：
[=ysum_{i=1}^nvarphi(frac{w}{gcd(w,i)})varphi(i)sum_{e|gcd(w,i)}varphi(e)
]
因为 (frac{w}{gcd(w,i)}) 和 (gcd(w,i)) 互质，所以他和 (e) 也互质，我们把他们结合：
[=ysum_{i=1}^nvarphi(i)sum_{e|gcd(w,i)}varphi(frac{w}{frac{gcd(w,i)}{e}})=ysum_{i=1}^nvarphi(i)sum_{e|gcd(w,i)}varphi(frac{w}{e})
]
可以先枚举 (e) 了，根据莫比乌斯反演的经验这样做一定有用：
[=ysum_{e|n}varphi(frac{w}{e})sum_{i=1}^{m/e}varphi(ie)=ysum_{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
",(ans+MOD)%MOD);
}