题目描述
求
[sum_{i=1}^nsum_{j=1}^msum_{k=1}^pgcd(icdot j,icdot k,jcdot k) imes gcd(i,j,k) imes left(frac{gcd(i,j)}{gcd(i,k) imes gcd(j,k)}+frac{gcd(i,k)}{gcd(i,j) imes gcd(j,k)}+frac{gcd(j,k)}{gcd(i,j) imes gcd(i,k)}
ight)
]
由于答案可能过大，输出答案对10^9+7109+7取模的值。
输入输出
第一行一个正整数T，为数据组数。
下面T行，每行3个整数，为n,m,p。
输出格式
共T行，每行一个整数，为答案。
输入输出
2
10 12 11
30 20 25
输出样例
25302
573830
Solution
挺巧妙的一个题。
注意到(gcd)的一个性质，我们只考虑一个质因子，设(i=p^x,j=p^y,k=p^z)，可以得到：
[gcd(i,j)=p^{min(x,y)},gcd(i,j,k)=p^{min(x,y,z)}
]
那么根据这个我们可以尝试着化简题目给出的式子，通分之后把分母提出来，和前面两项乘起来就是：
[res=dfrac{gcd(i,j,k) imes gcd(icdot j,jcdot k,kcdot i)}{gcd(i,j) imesgcd(j,k) imesgcd(i,k)}
]
由于这里只考虑一个质因子，我们可以两边取(log)，然后设(min(x,y,z)=x)，即(x)为最小值，这对答案是没有影响的，那么式子可以变成这样：
[log_p(res)=2x+min(y,y+z-x,z)-(2x+min(y,z))
]
可以发现(y+z-xgeqslant y,y+z-xgeqslant z)，所以可以得到：
[log_p(res)=2x+min(y,z)-(2x+min(y,z))
]
所以：
[log_p(res)=0,res=1
]
所以我们可以惊奇的发现，前面两项和分母约掉了，剩下的式子写出来就是：
[ans=sum_{i}sum_{j}sum_{k}gcd(i,j)^2+gcd(j,k)^2+gcd(i,k)^2
]
这个直接大力反演就好了。
#include<bits/stdc++.h>
using namespace std;
void read(int &x) {
x=0;int f=1;char ch=getchar();
for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}
void print(int x) {
if(x<0) putchar('-'),x=-x;
if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('
');}
const int maxn = 2e7+10;
const int mod = 1e9+7;
int pri[maxn/10],mu[maxn],f[maxn],vis[maxn],tot;
void sieve() {
mu[1]=f[1]=1;
for(int i=2;i<maxn;i++) {
if(!vis[i]) pri[++tot]=i,mu[i]=-1,f[i]=1ll*i*i%mod-1;
for(int v,j=1;j<=tot&&i*pri[j]<maxn;j++) {
vis[v=i*pri[j]]=1;
if(i%pri[j]==0) {f[v]=1ll*f[i]*pri[j]%mod*pri[j]%mod;break;}
f[v]=1ll*f[i]*f[pri[j]]%mod,mu[v]=-mu[i];
}
}
for(int i=1;i<maxn;i++) f[i]=(f[i]+f[i-1])%mod;
}
int calc(int n,int m) {
int T=1,res=0;
while(T<=min(n,m)) {
int pre=T;T=min(n/(n/T),m/(m/T));
res=(res+1ll*(f[T]-f[pre-1])*(n/T)%mod*(m/T)%mod)%mod;T++;
}return (res+mod)%mod;
}
int main() {
sieve();
int T,n,m,p;read(T);
while(T--) read(n),read(m),read(p),write(((1ll*calc(n,m)*p%mod+1ll*calc(n,p)*m%mod)%mod+1ll*calc(m,p)*n%mod)%mod);
return 0;
}