【LGP5176】公约数

题目

老年选手啥都不想推只能搞个杜教筛了

这个式子长得好吓人啊，所以我们唯一分解之后来考虑这道题

设\(i,j,k\)分别是\(p^a,p^b,p^c\)，至于到底谁是谁并不重要，我们不妨假设\(a\leq b\leq c\)

那么\(gcd(i\times j,j\times k,k\times 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}^n\sum_{j=1}^m\sum_{k=1}^p(i,j)^2+(j,k)^2+(k,i)^2 \]

显然可以拆成

\[p\times \sum_{i=1}^n\sum_{j=1}^m(i,j)^2+n\times \sum_{j=1}^m\sum_{k=1}^p(j,k)^2+m\times \sum_{i=1}^n\sum_{k=1}^p(k,i)^2 \]

套路反演我们可以得到我们要求的东西实际上是

\[\sum_{i=1}^n\left \lfloor \frac{n}{i} \right \rfloor\left \lfloor \frac{m}{i} \right \rfloor\sum_{d|i}\mu(\frac{i}{d})d^2 \]

后面的东西线筛也很好推，但是老年选手并不想动脑子了

发现后面的柿子是\(id^2\times \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\n",tmp%mod);
	}
}