[luogu5176] 公约数

Posted 2021-02-12 hbyer

题目描述

求
[ 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-xleqslant y,y+z-xleqslant 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;
}