Luogu5176 公约数 | Arahc's Home

题面

简化题意

设 $f(i,j,k) = \frac{\gcd(i,j)}{\gcd(i,k)\gcd(j,k)}$ ，求下面式子的值，对 $10^9+7$ 取模：

\sum_{i=1}^n \sum_{j=1}^m \sum_{k=1}^p \gcd(ij,ik,jk)\gcd(i,j,k)(f(i,j,k)+f(i,k,j)+f(j,k,i))

多测， $T\leq 10^3,n,m,p\leq 2\times 10^7$ 。

题解

考虑唯一分解定理，将 $\gcd(ij,jk,ik)$ 先干掉：将 $i,j,k$ 两两的公共质因子的并，去掉 $i,j,k$ 都有的质因子，剩下的部分就是 $ij,jk,ik$ 的质因子的交。也就是：

\gcd(ij,ik,jk)\gcd(i,j,k) = \gcd(i,j)\gcd(i,k)\gcd(i,k)

接下来考虑三个 $f$ 的和。设 $a=\gcd(i,j),b=\gcd(i,k),c=\gcd(j,k)$ ，则那三个和可以写成 $\frac{a}{bc}+\frac{b}{ac}+\frac{c}{ab}$ 。通分一下，代回原式：

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

设 $F(n,m,p) = p \sum_i^n \sum_j^m \gcd(i,j)^2$ 。显然只要能求出一个 $F$ 怎么算原问题就解决里。

考虑后面的 $\gcd(i,j)^2$ 很莫反：

\sum_{d=1}^{\min(n,m)} d^2 \sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor} \sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor} \sum_{k\mid \gcd(i,j)} \mu(k)

然后交换 $\sum$ ：

\sum_{d=1}^{\min(n,m)} d^2 \sum_{k=1}^{\left\lfloor\frac{\min(n,m)}{d}\right\rfloor} \mu(k) \left\lfloor\frac{n}{dk}\right\rfloor \left\lfloor\frac{m}{dk}\right\rfloor

然后换元，设 $t=dk$ ：

\sum_{p=1}^n \left\lfloor\frac{n}{p}\right\rfloor \left\lfloor\frac{m}{p}\right\rfloor \sum_{d\mid p} d^2\mu(\frac{p}{d})

整除分块，后面是一个 $\rm{id}^2\ast\mu$ 的形式。可以杜教筛。不过数据范围开得比较小，线性筛就能过了。

代码

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int max_n=20000007,N=20000000,mod=1000000007;
inline int read(){
    int x=0;bool w=0;char c=getchar();
    while(c<'0' || c>'9') w|=c=='-',c=getchar();
    while(c>='0' && c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return w?-x:x;
}
inline void write(int x){
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10^48);
}

bool isp[max_n];
int pr[max_n],cntp,mu[max_n],f[max_n];

inline void Getpr(int n){
    memset(isp,1,sizeof(isp)),isp[1]=0;
    mu[1]=f[1]=1;
    for(register int i=2;i<=n;++i){
        if(isp[i]){
            pr[++cntp]=i;
            mu[i]=-1,
            f[i]=(i*i-1)%mod;
        }
        for(register int j=1;i*pr[j]<=n && j<=cntp;++j){
            isp[i*pr[j]]=0;
            if(i%pr[j]==0){
                mu[i*pr[j]]=0,
                f[i*pr[j]]=f[i]*pr[j]*pr[j]%mod;
                break;
            }
            mu[i*pr[j]]=-mu[i],
            f[i*pr[j]]=f[i]*f[pr[j]]%mod;
        }
    }
    for(register int i=1;i<=n;++i)
        f[i]=(f[i-1]+f[i])%mod;
}

inline int F(int n,int m){
    int res=0,up=min(n,m);
    for(register int l=1,r;l<=up;l=r+1){
        r=min(n/(n/l),m/(m/l));
        res=(res+(n/l)*(m/l)%mod*(f[r]-f[l-1]+mod)%mod)%mod;
    }
    return res;
}

signed main(){
    Getpr(N);
    for(register int T=read();T>0;--T){
        int n=read(),m=read(),p=read();
        write((p*F(n,m)+m*F(n,p)+n*F(m,p))%mod),putchar('\n');
    }
    return 0;
}