GCD(i,j)求和（莫比乌斯函数）（CJOJ2512）

最新推荐文章于 2025-05-03 18:27:35 发布

JK Chen

最新推荐文章于 2025-05-03 18:27:35 发布

阅读量5.7k

点赞数 3

CC 4.0 BY-SA版权

分类专栏：数论/数学例题

本文链接：https://blog.youkuaiyun.com/jk_chen_acmer/article/details/82712276

数论/数学例题专栏收录该内容

70 篇文章

订阅专栏

本文介绍了一种解决双层求和中最大公约数(GCD)问题的有效算法。通过转换视角到因子角度并利用莫比乌斯反演，将问题简化为可快速求解的形式。适用于n,m≤1e7的大规模数据处理。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题意：

求 $∑i=1n∑j=1mgcd(i,j)\sum_{i=1}^n\sum_{j=1}^mgcd(i,j)$
n,m=1e7

$O (N)$ 版本：

换成因子角度，设： $g c d (i, j) = d, m i n (n, m) = N$ ，有： $Ans=∑d=1Nd∑i=1n∑j=1m[gcd(i,j)==d]=∑d=1Nd∑i=1nd∑j=1md[gcd(i,j)==1]Ans=\sum_{d=1}^Nd\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)==d]=\sum_{d=1}^Nd\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}[gcd(i,j)==1]$ 而 $g c d = = d$ 按照套路可以进行莫比乌斯的倍数反演： $g(d,n,m)=∑i=1n∑j=1m[gcd(i,j)==d]构造：f(d,n,m)=∑d∣Dg(D,n,m)=∑i=1n∑j=1m[d∣gcd(i,j)]=[nd][md]反演：g(d,n,m)=∑d∣Dmin(n,m)u(Dd)f(D)=∑d∣Dmin(n,m)u(Dd)[nD][mD]而我们所求为gcd==1，所以：g(1,n,m)=∑i=1Nu(i)[ni][mi]Ans=∑d=1Nd∗g(1,nd,md)Ans=∑d=1Nd∑i=1Ndu(i)[ni∗d][mi∗d]g(d,n,m)=\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)==d] \\ 构造：f(d,n,m)=\sum_{d|D}g(D,n,m)=\sum_{i=1}^n\sum_{j=1}^m[d|gcd(i,j)]=[\frac{n}{d}][\frac{m}{d}]\\反演：g(d,n,m)=\sum_{d|D}^{min(n,m)}u(\frac{D}{d})f(D)=\sum_{d|D}^{min(n,m)}u(\frac{D}{d})[\frac{n}{D}][\frac{m}{D}]\\而我们所求为gcd==1，所以：g(1,n,m)=\sum_{i=1}^Nu(i)[\frac{n}{i}][\frac{m}{i}]\\Ans=\sum_{d=1}^Nd*g(1,\frac{n}{d},\frac{m}{d})\\Ans=\sum_{d=1}^Nd\sum_{i=1}^\frac{N}{d}u(i)[\frac{n}{i*d}][\frac{m}{i*d}]$ 同样，到这里两个分块时间复杂度为O(N)

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const int  N = 1e7+9;
const LL mod = 998244353;
LL u[N];
LL pri[N>>1],now;
bool vis[N];
void init(){                            //预处理莫比乌斯函数前缀和
    vis[1]=1;u[1]=1;
    for(int i=2;i<N;i++){
        if(!vis[i])pri[++now]=i,u[i]=-1;
        for(int j=1;j<=now&&pri[j]*i<N;j++){
            vis[pri[j]*i]=1;
            u[pri[j]*i]=-u[i];
            if(i%pri[j]==0){u[pri[j]*i]=0;break;}
        }
    }
    for(int i=1;i<N;i++)u[i]=u[i]+u[i-1];
}

LL solve(LL n,LL m){
    LL ans=0;
    LL NN=min(n,m);
    for(LL l=1,r;l<=NN;l=r+1){          //[(n/d)/i][(m/d)/i] 分块
        r=min(n/(n/l),m/(m/l));
        LL re=(u[r]-u[l-1])*(n/l)%mod*(m/l)%mod ;
        ans=(ans+re)%mod;
    }
    return (ans+mod)%mod;
}

int main(){
    init();
    LL n,m;
    while(cin>>n>>m){
        LL NN=min(n,m);
        LL ans=0;
        for(LL l=1,r;l<=NN;l=r+1){      //N/d 分块
            r=NN/(NN/l);
            LL d=(l+r)*(r-l+1)/2%mod;
            LL re=solve(n/l,m/l);
            re=(re+mod)%mod;
            ans=(ans+d*re%mod)%mod;
        }
        printf("%lld\n",ans);
    }
}