洛谷P2522 - [HAOI2011]Problem b

本文深入探讨了莫比乌斯反演算法在解决特定数学问题中的应用,通过实例详细讲解了如何利用该算法高效计算两数最大公约数等于特定值的数对个数。文章提供了一种优化计算策略,复杂度为O(T√10^5),并附带源代码实现。

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

Portal

Description

进行\(T(T\leq10^5)\)次询问,每次给出\(x_1,x_2,y_1,y_2\)\(d\)(均不超过\(10^5\)),求\(\sum_{i=x_1}^{x_2} \sum_{j=y_1}^{y_2} [gcd(i,j)=d]\)

Solution

莫比乌斯反演入门题。
\(calc(n,m)\)表示\(i\in[1,n],j\in[1,m]\)\(gcd(i,j)=d\)的数对\((i,j)\)的个数。那么简单地进行容斥,可知\(ans=calc(x_2,y_2)-calc(x_1-1,y_2)-calc(x_2,y_1-1)+calc(x_1-1,x_2-1)\)
于是考虑如何计算\(calc(n,m)\)
\[ f(d) = \sum_{i=1}^n \sum_{j=1}^m [gcd(i,j)=d] \]

\[\begin{align*} F(x) &= \sum_{x|d} f(d) \\ &= \sum_{x|d} \sum_{i=1}^n \sum_{j=1}^m [gcd(i,j)=d] \\ &= \sum_{k=1}^{⌊\frac{n}{x}⌋} \sum_{i=1}^n \sum_{j=1}^m [gcd(i,j)=kx] \\ &= ⌊\frac{n}{x}⌋⌊\frac{m}{x}⌋ \end{align*}\] \(gcd(i,j)=kx \Leftrightarrow x|i\)\(x|j\),那满足条件的\((i,j)\)就有\(⌊\frac{n}{x}⌋⌊\frac{m}{x}⌋\)对。再进行莫比乌斯反演:
\[ f(x)= \sum_{x|d} \mu(\frac{d}{x}) F(d) = \sum_{x|d} \mu(\frac{d}{x})⌊\frac{n}{d}⌋⌊\frac{m}{d}⌋ = \sum_{k=1}^{⌊\frac{n}{x}⌋} \mu(k)⌊\frac{n}{kx}⌋⌊\frac{m}{kx}⌋ \]这个做法看起来是\(O(\dfrac{n}{x})\)的。不过由于\(⌊\dfrac{n}{i}⌋\)最多只有\(\sqrt n\)种取值,所以我们可以以\(O(\sqrt n)\)的复杂度进行计算。

i123456789101112131415
15/i1575332211111111

观察发现,一个取值为\(v\)的区间是以\(⌊\frac{n}{v}⌋\)结尾的,下一个区间是从\(⌊\frac{n}{v}⌋+1\)开始的,模拟这一性质去计算即可。若对于区间\(k\in[L,R]\)\(⌊\frac{n}{kx}⌋=v_1,⌊\frac{m}{kx}⌋=v_2\),那么该区间对答案的贡献为\(v_1v_2\sum_{k=L}^R \mu(k)\),预处理出\(\mu(x)\)的前缀和即可。

时间复杂度\(O(T\sqrt {10^5})\)

Code

//[HAOI2011]Problem b
#include <algorithm>
#include <cstdio>
using std::min; using std::swap;
typedef long long lint;
inline char gc()
{
    static char now[1<<16],*s,*t;
    if(s==t) {t=(s=now)+fread(now,1,1<<16,stdin); if(s==t) return EOF;}
    return *s++;
}
inline int read()
{
    int x=0; char ch=gc();
    while(ch<'0'||'9'<ch) ch=gc();
    while('0'<=ch&&ch<='9') x=x*10+ch-'0',ch=gc();
    return x;
}
const int N=5e4+10;
int mu[N],pre[N];
int cntP,pr[N]; bool notP[N];
void getMu(int n)
{
    mu[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!notP[i]) pr[++cntP]=i,mu[i]=-1;
        for(int j=1;j<=cntP;j++)
        {
            if((lint)i*pr[j]>n) break;
            int x=i*pr[j]; notP[x]=true;
            if(i%pr[j]) mu[x]=-mu[i]; else {mu[x]=0; break;}
        }
    }
    for(int i=1;i<=n;i++) pre[i]=pre[i-1]+mu[i];
}
int k;
lint calc(int x,int y)
{
    x/=k,y/=k; if(x>y) swap(x,y);
    lint res=0;
    for(int L=1,R;L<=x;L=R+1)
    {
        int v1=x/L,v2=y/L; R=min(x/v1,y/v2);
        res+=1LL*(pre[R]-pre[L-1])*v1*v2;
    }
    return res;
}
int main()
{
    getMu(5e4);
    int Q=read();
    while(Q--)
    {
        int fr1=read(),to1=read(),fr2=read(),to2=read(); k=read();
        printf("%lld\n",calc(to1,to2)-calc(fr1-1,to2)-calc(to1,fr2-1)+calc(fr1-1,fr2-1));
    }
    return 0;
}

P.S.

同样的题洛谷P2257

转载于:https://www.cnblogs.com/VisJiao/p/LgP2522.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值