题链：
http://www.lydsy.com/JudgeOnline/problem.php?id=3512
题解：
$$求ANS=sum_{i=1}^{N}sum_{j=1}^{M}phi(ij)quad Nleq 10^5;Mleq 10^9$$
杜教筛
因为N比较小，所以从这里入手：
设$sum(n,M)=sum_{i=1}^{M}phi(ni)$
则答案为$ANS=sum_{n=1}^{N}sum(n,M)$
考虑如何求$sum(n,M)$
首先按照唯一分解定理，$n={p_1}^{a_1} imes {p_2}^{a_2} imescdots imes {p_k}^{a_k}$
另$n'={p_1}^{1} imes {p_2}^{1} imescdots imes {p_k}^{1}，P=frac{n}{n'}$
由$phi$的定义可得：$$phi(n)=Pphi(frac{n}{P})=Pphi(n')$$
所以$$sum(n,M)=P imes sum(n',M)$$
（现在n'没有平方质因子，即$|mu(n')|=1$）
第一阶段已经结束，我们看看$sum(n',M)$又该怎么求？
任取一个n'的质因子，我们这里取最小的那个，用miniP表示：
显然$gcd(miniP,frac{n'}{miniP})=1，即miniP和frac{n'}{miniP}互质$
从$sum(n',M)=sum_{i=1}^{M}phi(n'i)$这个定义入手，接下来分两种情况：
（一）、miniP不是i的因子
那么有$$phi(n'i)=phi(miniP)phi(frac{n'}{miniP}i)=(miniP-1)phi(frac{n'}{miniP}i)$$
先假设$sum(n',M)=sum_{i=1}^{M}phi(n'i)$中枚举的i都与miniP互质，那么得到：
$$egin{aligned}
sum(n',M)&=sum_{i=1}^{M}phi(n'i)\
&=sum_{i=1}^{M}(miniP-1)phi({frac{n'}{minip}i})\
&=(miniP-1)sum(frac{n'}{minip},M)\
end{aligned}$$
（二）、miniP是i的因子
那么有$$phi(n'i)=miniPphi(frac{n'}{miniP}i)$$
显然对于上面的$sum(n',M)=(miniP-1)sum(frac{n'}{minip},M)$而言，每当枚举到一个i与miniP不互质时，
就会少加一个$phi(frac{n'}{miniP}i)$
现在我们希望把这些漏掉的加回来，设漏加的总和为w，则：
$$egin{aligned}
w&=sum_{miniP|i}phi(frac{n'}{miniP}i)\
&=sum_{i=1}^{lfloor frac{M}{miniP}
floor}phi(frac{n'}{miniP}(i imes miniP))\
&=sum_{i=1}^{lfloor frac{M}{miniP}
floor}phi(n'i)\
&=sum(n',lfloor frac{M}{miniP}
floor)
end{aligned}$$
所以，综上两种情况，我们得到：
$$sum(n',M)=(miniP-1)sum(frac{n'}{minip},M)+sum(n',lfloor frac{M}{miniP} floor)$$
差不多就这样了，我们枚举每个n，
分别求出：
$sum(n,M)=P imes sum(n',M)=P imes((miniP-1)sum(frac{n'}{minip},M)+sum(n',lfloor frac{M}{miniP} floor))$
对于每个sum()，进行递归求解，
到了递归的最底层时，
若n==1，那么$sum(1,M)=sum_{i=1}^{M}phi(i)$，用杜教筛求解。
若m==1，那么$sum(n,1)=phi(n)，直接返回即可$
另外要先预处理出前$M^{frac{2}{3}}个phi()$的前缀和，便于降低杜教筛的复杂度。
代码：
#include<bits/stdc++.h>
#define MAXN 100050
#define DJM 1000000
#define mod 1000000007
using namespace std;
struct Hash_Table{
#define Hmod 1425367
int org[DJM],val[DJM],nxt[DJM],head[Hmod],hnt;
Hash_Table(){hnt=2;}
void Push(int x,int v){
static int u; u=x%Hmod;
org[hnt]=x; val[hnt]=v; nxt[hnt]=head[u]; head[u]=hnt++;
}
int Find(int x){
static int u; u=x%Hmod;
for(int i=head[u];i;i=nxt[i])
if(org[i]==x) return val[i];
return -1;
}
}H;
int P[DJM+50],miniP[DJM+50],phi[DJM+50],sphi[DJM+50];
int sum[MAXN];
int N,M;
void Sieve(){
static bool np[DJM+50];
static int prime[DJM+50],pnt;
phi[1]=1; P[1]=1; miniP[1]=1;
for(int i=2;i<=DJM;i++){
if(!np[i]) prime[++pnt]=i,phi[i]=i-1,P[i]=1,miniP[i]=i;
for(int j=1;j<=pnt&&i<=DJM/prime[j];j++){
np[i*prime[j]]=1; miniP[i*prime[j]]=prime[j];
if(i%prime[j]){
phi[i*prime[j]]=phi[i]*phi[prime[j]];
P[i*prime[j]]=P[i];
}else{
phi[i*prime[j]]=phi[i]*prime[j];
P[i*prime[j]]=P[i]*prime[j];
break;
}
}
}
//for(int i=1;i<=DJM;i++) printf(":P=%d miniP=%d phi=%d
",P[i],miniP[i],phi[i]);
for(int i=1;i<=DJM;i++) sphi[i]=(1ll*phi[i]+sphi[i-1])%mod;
}
int DJ_Sieve(int m){
if(m<=DJM) return sphi[m];
if(H.Find(m)!=-1) return H.Find(m);
int now=(1ll*m*(1+m)/2)%mod;
for(int i=2,last;i<=m;i=last+1){
last=m/(m/i);
now=(1ll*now-1ll*(last-i+1)*DJ_Sieve(m/i)%mod+mod)%mod;
}
H.Push(m,now);
return now;
}
int S(int n,int m){
if(m==0) return 0;
if(m==1) return phi[n];
if(n==1) return DJ_Sieve(m);
if(m==M&&sum[n]) return sum[n];
int now=(1ll*(miniP[n]-1)*S(n/miniP[n],m)+S(n,m/miniP[n]))%mod;
if(m==M) sum[n]=now;
if(now<0) printf("S %d %d
",n,m);
return now;
}
int main(){
Sieve();
scanf("%d%d",&N,&M);
int ans=0;
for(int n=1;n<=N;n++)
ans=(1ll*ans+1ll*P[n]*S(n/P[n],M))%mod;
printf("%d
",ans);
return 0;
}