题链：
题解：
$$求ANS=\sum_{i=1}^{N}\sum_{j=1}^{M}\phi(ij)\quad N\leq 10^5\;M\leq 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}\times {p_2}^{a_2}\times\cdots\times {p_k}^{a_k}$
另$n'={p_1}^{1}\times {p_2}^{1}\times\cdots\times {p_k}^{1}，P=\frac{n}{n'}$
由$\phi$的定义可得：$$\phi(n)=P\phi(\frac{n}{P})=P\phi(n')$$
所以$$sum(n,M)=P\times 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互质，那么得到：
$$\begin{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)=miniP\phi(\frac{n'}{miniP}i)$$
显然对于上面的$sum(n',M)=(miniP-1)sum(\frac{n'}{minip},M)$而言，每当枚举到一个i与miniP不互质时，
就会少加一个$\phi(\frac{n'}{miniP}i)$
现在我们希望把这些漏掉的加回来，设漏加的总和为w，则：
$$\begin{aligned}w&=\sum_{miniP|i}\phi(\frac{n'}{miniP}i)\\&=\sum_{i=1}^{\lfloor \frac{M}{miniP} \rfloor}\phi(\frac{n'}{miniP}(i\times miniP))\\&=\sum_{i=1}^{\lfloor \frac{M}{miniP} \rfloor}\phi(n'i)\\&=sum(n',\lfloor \frac{M}{miniP} \rfloor)\end{aligned}$$
所以，综上两种情况，我们得到：
$$sum(n',M)=(miniP-1)sum(\frac{n'}{minip},M)+sum(n',\lfloor \frac{M}{miniP} \rfloor)$$
差不多就这样了，我们枚举每个n，
分别求出：
$sum(n,M)=P\times sum(n',M)=P\times((miniP-1)sum(\frac{n'}{minip},M)+sum(n',\lfloor \frac{M}{miniP} \rfloor))$
对于每个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#define MAXN 100050#define DJM 1000000#define mod 1000000007using 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\n",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",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\n",ans); return 0;}