【BZOJ3512】DZY Loves Math IV

题目大意

给出 n , m n,m n,m，求 ∑ i = 1 n ∑ j = 1 m ϕ ( i j ) \sum\limits_{i=1}^n\sum\limits_{j=1}^m\phi(ij) i=1∑nj=1∑mϕ(ij)

1 ≤ n ≤ 1 0 5 , 1 ≤ m ≤ 1 0 9 1\leq n\leq 10^5,1\leq m\leq 10^9 1≤n≤105,1≤m≤109，答案模 1 0 9 + 7 10^9+7 109+7。

题解

前置知识：杜教筛

因为 n n n比较小，所以我们可以直接枚举 n n n。

设 f ( n , m ) = ∑ i = 1 m ϕ ( n i ) f(n,m)=\sum\limits_{i=1}^m\phi(ni) f(n,m)=i=1∑mϕ(ni)

将 n n n分为若干个质数的乘积， n = p 1 a 1 p 2 a 2 ⋯ p k a k n=p_1^{a_1}p_2^{a_2}\cdots p_k^{a_k} n=p1a1p2a2⋯pkak，令 p = p 1 a 1 − 1 p 2 a 2 − 1 ⋯ p k a k − 1 p=p_1^{a_1-1}p_2^{a_2-1}\cdots p_k^{a_k-1} p=p1a1−1p2a2−1⋯pkak−1， q = n / p q=n/p q=n/p

∑ i = 1 m ϕ ( n i ) = ∑ i = 1 m p ϕ ( q i ) = ∑ i = 1 m p ϕ ( q gcd ⁡ ( q , i ) × i ) × gcd ⁡ ( q , i ) \sum\limits_{i=1}^m\phi(ni)=\sum\limits_{i=1}^mp\phi(qi)=\sum\limits_{i=1}^mp\phi(\dfrac{q}{\gcd(q,i)}\times i)\times \gcd(q,i) i=1∑mϕ(ni)=i=1∑mpϕ(qi)=i=1∑mpϕ(gcd(q,i)q×i)×gcd(q,i)

因为 ϕ ( i ) \phi(i) ϕ(i)为积性函数，所以 ϕ ( q gcd ⁡ ( q , i ) × i ) = ϕ ( q gcd ⁡ ( q , i ) ) × ϕ ( i ) \phi(\dfrac{q}{\gcd(q,i)}\times i)=\phi(\dfrac{q}{\gcd(q,i)})\times\phi(i) ϕ(gcd(q,i)q×i)=ϕ(gcd(q,i)q)×ϕ(i)，原式变为

p ∑ i = 1 m ϕ ( q gcd ⁡ ( q , i ) ) × ϕ ( i ) × gcd ⁡ ( q , i ) p\sum\limits_{i=1}^m\phi(\dfrac{q}{\gcd(q,i)})\times\phi(i)\times \gcd(q,i) pi=1∑mϕ(gcd(q,i)q)×ϕ(i)×gcd(q,i)

由欧拉函数的性质， n = ∑ d ∣ n ϕ ( n d ) n=\sum\limits_{d|n}\phi(\dfrac nd) n=d∣n∑ϕ(dn)，可得

p ∑ i = 1 m ϕ ( i ) ϕ ( q gcd ⁡ ( q , i ) ) ∑ j ∣ i , j ∣ q ϕ ( gcd ⁡ ( q , i ) j ) p\sum\limits_{i=1}^m\phi(i)\phi(\dfrac{q}{\gcd(q,i)})\sum\limits_{j|i,j|q}\phi(\dfrac{\gcd(q,i)}{j}) pi=1∑mϕ(i)ϕ(gcd(q,i)q)j∣i,j∣q∑ϕ(jgcd(q,i))

因为 q gcd ⁡ ( q , i ) \dfrac{q}{\gcd(q,i)} gcd(q,i)q与 gcd ⁡ ( q , i ) \gcd(q,i) gcd(q,i)互质（ q = p 1 p 2 ⋯ p k q=p_1p_2\cdots p_k q=p1p2⋯pk，同一个质数要不不出现，要不出现在 q gcd ⁡ ( q , i ) \dfrac{q}{\gcd(q,i)} gcd(q,i)q或 gcd ⁡ ( q , i ) \gcd(q,i) gcd(q,i)，不可能两个都出现）而 ϕ \phi ϕ为积性函数，所以 ϕ ( q gcd ⁡ ( q , i ) ) × ϕ ( gcd ⁡ ( q , i ) j ) = ϕ ( q j ) \phi(\dfrac{q}{\gcd(q,i)})\times \phi(\dfrac{\gcd(q,i)}{j})=\phi(\dfrac qj) ϕ(gcd(q,i)q)×ϕ(jgcd(q,i))=ϕ(jq)，可得

p ∑ i = 1 m ϕ ( i ) ∑ j ∣ i , j ∣ q ϕ ( q j ) p\sum\limits_{i=1}^m\phi(i)\sum\limits_{j|i,j|q}\phi(\dfrac qj) pi=1∑mϕ(i)j∣i,j∣q∑ϕ(jq)

枚举 j j j得

p ∑ j ∣ q ϕ ( q j ) ∑ i = 1 ⌊ m j ⌋ ϕ ( i j ) p\sum\limits_{j|q}\phi(\dfrac qj)\sum\limits_{i=1}^{\lfloor\frac mj\rfloor}\phi(ij) pj∣q∑ϕ(jq)i=1∑⌊jm⌋ϕ(ij)

因为 f ( n , m ) = ∑ i = 1 m ϕ ( n i ) f(n,m)=\sum\limits_{i=1}^m\phi(ni) f(n,m)=i=1∑mϕ(ni)，所以

p ∑ j ∣ q ϕ ( q j ) × f ( j , ⌊ m j ⌋ ) p\sum\limits_{j|q}\phi(\dfrac qj)\times f(j,\lfloor\dfrac mj\rfloor) pj∣q∑ϕ(jq)×f(j,⌊jm⌋)

于是我们就得到了递推式

f ( n , m ) = p ∑ j ∣ q ϕ ( q j ) × f ( j , ⌊ m j ⌋ ) f(n,m)=p\sum\limits_{j|q}\phi(\dfrac qj)\times f(j,\lfloor\dfrac mj\rfloor) f(n,m)=pj∣q∑ϕ(jq)×f(j,⌊jm⌋)

预处理出前 n 2 3 n^{\frac 23} n32个 ϕ ( i ) \phi(i) ϕ(i)并求前缀和，用杜教筛 ϕ ( i ) \phi(i) ϕ(i)的前缀和，利用 ϕ \phi ϕ的前缀和作差得出每个 ϕ ( q j ) \phi(\dfrac qj) ϕ(jq)，然后就可以用杜教筛求 f ( n , m ) f(n,m) f(n,m)了。

注意 f ( n , m ) f(n,m) f(n,m)的一些特殊值。

当 n = 1 n=1 n=1时， f ( n , m ) = ∑ i = 1 m ϕ ( i ) f(n,m)=\sum\limits_{i=1}^m\phi(i) f(n,m)=i=1∑mϕ(i)
当 m = 1 m=1 m=1时， f ( n , m ) = ϕ ( n ) f(n,m)=\phi(n) f(n,m)=ϕ(n)

最终的答案为 ∑ i = 1 n f ( i , m ) \sum\limits_{i=1}^nf(i,m) i=1∑nf(i,m)。

code

#include<bits/stdc++.h>
using namespace std;
const int N=2000000;
const long long mod=1000000007;
int z[N+5],p[N+5],mn[N+5];
long long ans,ph[N+5],s[N+5];
map<long long,long long>sph,sum[100005];
void init(){ph[1]=1;for(int i=2;i<=N;i++){if(!z[i]){p[++p[0]]=i;ph[i]=i-1;mn[i]=i;}for(int j=1;j<=p[0]&&i*p[j]<=N;j++){z[i*p[j]]=1;mn[i*p[j]]=p[j];if(i%p[j]==0){ph[i*p[j]]=ph[i]*p[j]%mod;break;}ph[i*p[j]]=ph[i]*(p[j]-1)%mod;}}for(int i=1;i<=N;i++){s[i]=(s[i-1]+ph[i])%mod;}
}
long long gt(int n){if(n<=N) return s[n];if(sph.count(n)) return sph[n];long long re=0;for(int l=2,r;l<=n;l=r+1){r=min(n,n/(n/l));re=(re+1ll*(r-l+1)*gt(n/l)%mod)%mod;}re=(1ll*n*(n+1)/2-re+mod)%mod;return sph[n]=re;
}
long long djs(int n,int m){if(!m) return 0;if(n==1) return gt(m);if(m==1) return ph[n];if(sum[n].count(m)) return sum[n][m];long long re=0;vector<int>v;int p=1,q=1,t=n,l=0;while(t>1){int x=mn[t];q*=x;t/=x;v.push_back(x);++l;while(t%x==0){p*=x;t/=x;}}for(int i=0;i<1<<l;i++){int vt=1;for(int j=0;j<l;j++)if(i&(1<<j)) vt*=v[j];re=(re+ph[q/vt]*djs(vt,m/vt)%mod)%mod;}return sum[n][m]=re*p%mod;
}
int main()
{int n,m;init();scanf("%d%d",&n,&m);for(int i=1;i<=n;i++){ans=(ans+djs(i,m))%mod;}printf("%lld",ans);return 0;
}