老年选手啥都不想推只能搞个杜教筛了
这个式子长得好吓人啊，所以我们唯一分解之后来考虑这道题
设(i,j,k)分别是(p^a,p^b,p^c)，至于到底谁是谁并不重要，我们不妨假设(aleq bleq c)
那么(gcd(i imes j,j imes k,k imes i)=min{p^{a+b},p^{a+c},p^{b+c}}=p^{a+b})
(gcd(i,j,k)=min{p^a,p^b,p^c}=p^a)
也就是前面两个柿子乘起来是(p^{2a+b})
我们把后面的柿子分母通分
[frac{(i,j)^2+(j,k)^2+(k,i)^2}{(i,j)(j,k)(k,i)}
]
我们发现分母上还是等于(p^{2a+b})，因为(p^a)跟另外两个组合得到的(gcd)都是(p^a)，(p^b)和(p^{c})的(gcd)是(p^b)，所以分母上是(p^{2a+b})
和外面一约分，没了
显然唯一分解之后各个质数次幂是相互独立的，于是我们现在可以得出结论，我们要求的就是
[sum_{i=1}^nsum_{j=1}^msum_{k=1}^p(i,j)^2+(j,k)^2+(k,i)^2
]
显然可以拆成
[p imes sum_{i=1}^nsum_{j=1}^m(i,j)^2+n imes sum_{j=1}^msum_{k=1}^p(j,k)^2+m imes sum_{i=1}^nsum_{k=1}^p(k,i)^2
]
套路反演我们可以得到我们要求的东西实际上是
[sum_{i=1}^nleft lfloor frac{n}{i}
ight
floorleft lfloor frac{m}{i}
ight
floorsum_{d|i}mu(frac{i}{d})d^2
]
后面的东西线筛也很好推，但是老年选手并不想动脑子了
发现后面的柿子是(id^2 imes mu)，我们直接杜教筛卷上(1)就变成(id^2)了，(id^2)的前缀和自然是(frac{n(n+1)(2n+1)}{6})了
这样是(O(n^{frac{3}{4}}))的，但是我们暴力调合级数处理一些，就能跑的很快了
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define re register
#define LL long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read() {
char c=getchar();int x=0;while(c<'0'||c>'9') c=getchar();
while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();return x;
}
int T;
const int I6=166666668;
const int mod=1e9+7;
const int maxn=1e6+1;
int p[maxn>>1],is[maxn],mu[maxn];
int f[20000005],vis[20000005];
inline int calc(int x) {return 1ll*x*(x+1)%mod*(x+x+1)%mod*I6%mod;}
inline int solve(int x) {
if(vis[x]) return f[x];
vis[x]=1;int ans=calc(x);
for(re int l=2,r;l<=x;l=r+1) {
r=x/(x/l);
ans=(ans-1ll*solve(x/l)*(r-l+1)%mod+mod)%mod;
}
return f[x]=ans;
}
inline int work(int n,int m) {
if(n>m) std::swap(n,m);
int ans=0;
for(re int l=1,r;l<=n;l=r+1) {
r=min(n/(n/l),m/(m/l));
ans=(ans+1ll*(n/l)*(m/l)%mod*(solve(r)-solve(l-1)+mod)%mod)%mod;
}
return ans;
}
int main() {
T=read();is[1]=mu[1]=1;
for(re int i=2;i<maxn;i++) {
if(!is[i]) p[++p[0]]=i,mu[i]=-1;
for(re int j=1;j<=p[0]&&p[j]*i<maxn;j++) {
is[p[j]*i]=1;if(i%p[j]==0) break;mu[p[j]*i]=-1*mu[i];
}
}
for(re int i=1;i<maxn;i++)
for(re int j=i;j<maxn;j+=i) {
if(!mu[j/i]) continue;
if(mu[j/i]<0) f[j]=(f[j]-1ll*i*i%mod+mod)%mod;
else f[j]=(f[j]+1ll*i*i)%mod;
}
for(re int i=1;i<maxn;i++)
vis[i]=1,f[i]=(f[i]+f[i-1])%mod;
while(T--) {
int n=read(),m=read(),p=read();
int tmp=1ll*work(n,m)*p%mod+1ll*work(m,p)*n%mod;tmp%=mod;
tmp=(tmp+1ll*work(n,p)*m%mod)%mod;
printf("%d
",tmp%mod);
}
}