题目描述
\[\sum_{i=1}^n \sum_{j=1}^n \gcd(i,j)^2
\]
具体思路
solution 1
显然可以每次枚举 \(\gcd(i,j)\) 的取值。
\[\sum_{k=1}^n k^2 \sum_{i=1}^n \sum_{j=1}^n [\gcd(i,j)=k]
\]
令 \(i=\lfloor \frac{i}{k} \rfloor\)，\(j=\lfloor \frac{j}{k} \rfloor\)。
\[\sum_{k=1}^n k^2 \sum_{i=1}^{\lfloor \frac{n}{k} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{k} \rfloor} [\gcd(i,j)=1]
\]
\[\sum_{k=1}^n k^2 \sum_{i=1}^{\lfloor \frac{n}{k} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{k} \rfloor} \varepsilon(\gcd(i,j))
\]
然后进行莫比乌斯反演。
\[\sum_{k=1}^n k^2 \sum_{i=1}^{\lfloor \frac{n}{k} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{k} \rfloor} \sum_{d \mid \gcd(i,j)} \mu(d)
\]
\[\sum_{k=1}^n k^2 \sum_{i=1}^{\lfloor \frac{n}{k} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{k} \rfloor} \sum_{d \mid i,d \mid j} \mu(d)
\]
交换求和顺序，由于先枚举 \(i,j\) 再枚举 \(d\)，此时 \(d\) 是 \(i,j\) 的约数，因此反过来先枚举 \(d\) 再枚举 \(i,j\)，此时 \(i,j\) 就是 \(d\) 的倍数。而 \(1 \sim \lfloor \frac{n}{k} \rfloor\) 里面 \(d\) 的倍数有 \(\lfloor \frac{n}{kd} \rfloor\) 个。
\[\sum_{k=1}^n k^2 \sum_{d=1}^{\lfloor \frac{n}{k} \rfloor} \mu(d) \sum_{i=1}^{\lfloor \frac{n}{kd} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{kd} \rfloor} 1
\]
\[\sum_{k=1}^n k^2 \sum_{d=1}^{\lfloor \frac{n}{k} \rfloor} \mu(d) \lfloor \frac{n}{kd} \rfloor \lfloor \frac{n}{kd} \rfloor
\]
令：
\[f(n)=\sum_{d=1}^n \mu(d) \lfloor \frac{n}{d} \rfloor \lfloor \frac{n}{d} \rfloor
\]
有：
\[\sum_{k=1}^n k^2 f(\lfloor \frac{n}{k} \rfloor)
\]
显然 \(f(n)\) 可以数论分块解决，而最后的式子也可以数论分块解决，即数论分块套数论分块。
时间复杂度
\[O(\sum_{i=1}^{\sqrt n} \sqrt \frac{n}{i})=O(\sqrt n\int_0^{\sqrt n} i^{-\frac{1}{2}})=O(n^{\frac{3}{4}})
\]
这个复杂度跑 \(1e7\) 还是挺轻松的。
Code
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=1e7;
const int mod=1e9+7;
int v[N+5],prime[N+5],pr;
int mu[N+5];LL sum[N+5];
LL qpow(LL a,LL b){
LL ans=1%mod;a%=mod;
while(b){
if(b&1)ans=ans*a%mod;
a=a*a%mod;b>>=1;
}
return ans;
}
void init(int n){
memset(v,0,sizeof(v));pr=0;
mu[1]=1;
for(int i=2;i<=n;i++){
if(!v[i]){
prime[++pr]=i;
mu[i]=-1;
}
for(int j=1;j<=pr&&i*prime[j]<=n;j++){
v[i*prime[j]]=1;
if(i%prime[j]==0){
mu[i*prime[j]]=0;
break;
}
mu[i*prime[j]]=-mu[i];
}
}
for(int i=1;i<=n;i++){
sum[i]=(sum[i-1]+mu[i])%mod;
}
}
LL block(LL n){
LL ans=0;
for(LL l=1,r=0;l<=n;l=r+1){
r=n/(n/l);
ans=(ans+(sum[r]-sum[l-1])%mod*(n/l)%mod*(n/l)%mod+mod)%mod;
}
return ans;
}
LL f(LL n){
return n*(n+1)%mod*(2*n+1)%mod*qpow(6,mod-2);
}
LL calc(LL n){
LL ans=0;
for(LL l=1,r=0;l<=n;l=r+1){
r=n/(n/l);
ans=(ans+(f(r)-f(l-1))%mod*block(n/l)%mod+mod)%mod;
}
return ans;
}
int main(){
LL n;scanf("%lld",&n);
init(n);
printf("%lld\n",calc(n));
return 0;
}