[luogu5176] 公约数

题目描述

求

\[\sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p\gcd(i\cdot j,i\cdot k,j\cdot k)\times \gcd(i,j,k)\times \left(\frac{\gcd(i,j)}{\gcd(i,k)\times \gcd(j,k)}+\frac{\gcd(i,k)}{\gcd(i,j)\times \gcd(j,k)}+\frac{\gcd(j,k)}{\gcd(i,j)\times \gcd(i,k)}\right) \]

由于答案可能过大，输出答案对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)\times \gcd(i\cdot j,j\cdot k,k\cdot i)}{\gcd(i,j)\times\gcd(j,k)\times\gcd(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-x\geqslant y,y+z-x\geqslant 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('\n');}

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;
}