洛谷P2522 容斥原理+莫比乌斯反演

最新推荐文章于 2025-09-07 14:58:49 发布

原创最新推荐文章于 2025-09-07 14:58:49 发布 · 177 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#算法

这篇博客介绍了如何运用莫比乌斯反演解决涉及gcd(greatest common divisor)的计数问题，通过两个不同的方法，一是采用容斥原理结合莫比乌斯函数，二是直接在原始公式上进行操作，来高效地计算给定区间内满足gcd条件的整数对数量。这两种方法都涉及到数论分块和线性筛预处理，对于大规模数据具有较好的效率。

题意：

给出 $a, b, c, d, k$ ，计算
$\sum_{i=a}^{b}\sum_{j=c}^{d}[gcd(i,j)=k]$

Solution：

法一：

设
$f(n,m)=\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=k]$
容斥原理，原题即求
$f (b, d) - f (a - 1, d) - f (b, c - 1) + f (a - 1, c - 1)$
现在只需要求 $f (n, m)$ ，由于出现布尔表达式，考虑莫比乌斯反演，莫比乌斯反演需要 $[. . . = 1]$ 形式，于是改写成
$\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}[gcd(i,j)=1]$
莫比乌斯反演，即求
$\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}\sum_{d|gcd(i,j)}\mu(d)$
更换枚举项为因子 $d$
$\sum_{d=1}^{min\{\lfloor\frac{n}{k}\rfloor,\lfloor\frac{m}{k}\rfloor\}}\mu(d)\sum_{i=1}^{\lfloor\frac{n}{dk}\rfloor}\sum_{i=1}^{\lfloor\frac{m}{dk}\rfloor}1$
由于 $n, m$ 轮换，不妨设 $n≤mn\leq m$ ，即
$\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\lfloor\frac{n}{dk}\rfloor\lfloor\frac{m}{dk}\rfloor$
整除分块即可快速求得此式， $μ(d)\mu(d)$ 部分只需要线性筛后处理前缀和，需要注意的是，此处需要求单个 $μ(i)\mu(i)$ 而不是 $∑μ(i)\sum \mu(i)$ ，所以 $i%prime[j]=0i\%prime[j]=0$ 时结果是0

// #include<bits/stdc++.h>
#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<vector>
#include<bitset>
#include<map>
using namespace std;

using ll=long long;
const int N=5e4+5,inf=0x3fffffff;
const long long INF=0x3fffffffffffffff,mod=20101009;

bool nt[N];
int prime[N],cnt,mu[N];

void make_prime()
{
    mu[1]=1;
    for(int i=2;i<=50000;i++)
    {
        if(!nt[i]) prime[++cnt]=i,mu[i]=-1;
        for(int j=1;j<=cnt&&i*prime[j]<=50000;j++)
        {
            nt[i*prime[j]]=true;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            else mu[i*prime[j]]=mu[i]*mu[prime[j]];
        }
    }
    for(int i=1;i<=50000;i++) mu[i]+=mu[i-1];
}

ll solve(ll n,ll m,ll x)
{
    ll ret=0,l=1,r;
    n/=x; m/=x;
    while(l<=min(n,m))
    {
        r=min(min(n,m),min(n/(n/l),m/(m/l)));
        ret+=(n/l)*(m/l)*(mu[r]-mu[l-1]);
        l=r+1;
    }
    return ret;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    make_prime();
    int t; cin>>t;
    while(t--)
    {
        ll a,b,c,d,k; cin>>a>>b>>c>>d>>k;
        cout<<solve(b,d,k)-solve(a-1,d,k)-solve(b,c-1,k)+solve(a-1,c-1,k)<<'\n';
    }
    return 0;
}

法二：

不用容斥，直接在原式上操作
$\sum_{i=a}^{b}\sum_{j=c}^{d}[gcd(i,j)=k]\Rightarrow \sum_{i=\lceil\frac{a}{k}\rceil}^{\lfloor\frac{b}{k}\rfloor}\sum_{i=\lceil\frac{c}{k}\rceil}^{\lfloor\frac{d}{k}\rfloor}[gcd(i,j)=1]=\sum_{i=\lceil\frac{a}{k}\rceil}^{\lfloor\frac{b}{k}\rfloor}\sum_{i=\lceil\frac{c}{k}\rceil}^{\lfloor\frac{d}{k}\rfloor}\sum_{d|gcd(i,j)}\mu(d)$
注意上式的求和开始的下标是上取整

优先枚举因子 $d$ ，假设 $b≤db\leq d$ ， $f (x, l, r)$ 为 $x$ 在区间 $[l, r]$ 有多少个倍数
$\sum_{d=1}^{\lfloor\frac{b}{k}\rfloor}\mu(d)f(d,\lceil\frac{a}{k}\rceil,\lfloor\frac{b}{k}\rfloor)f(d,\lceil\frac{c}{k}\rceil,\lfloor\frac{d}{k}\rfloor)$
对于 $f (x, l, r)$ ，一个形式是
$f(x,l,r)=\lfloor\frac{r}{k}\rfloor-\lceil\frac{l}{k}\rceil+1$
这有上取整，是很不好处理的，下取整只能求得小于等于某个值的倍数有多少个，而 $[l, r]$ 的倍数个数可以容斥得到为小于等于 $r$ 的倍数的个数-小于等于 $l - 1$ 的倍数的个数，于是也可以写作如下形式
$f(x,l,r)=\lfloor\frac{r}{k}\rfloor-\lfloor\frac{l-1}{k}\rfloor$
代入原式后数论分块就可以求得答案，需要注意特判分母为0，不然会RE。这样的写法不需要容斥，常数小，更大的好处是容斥的时间复杂度是 $2^n$ ， $n$ 是容斥元素个数，当 $n$ 很大的时候是不能容斥的，需要这样子计算

// #include<bits/stdc++.h>
#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<vector>
#include<bitset>
#include<map>
using namespace std;

using ll=long long;
const int N=5e4+5,inf=0x3fffffff;
const long long INF=0x3fffffffffffffff,mod=20101009;

ll ceil(ll x,ll y){
    return x%y?x/y+1:x/y;
}

bool nt[N];
int prime[N],cnt,mu[N];

void make_prime()
{
    mu[1]=1;
    for(int i=2;i<=50000;i++)
    {
        if(!nt[i]) prime[++cnt]=i,mu[i]=-1;
        for(int j=1;j<=cnt&&i*prime[j]<=50000;j++)
        {
            nt[i*prime[j]]=true;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            else mu[i*prime[j]]=mu[i]*mu[prime[j]];
        }
    }
    for(int i=1;i<=50000;i++) mu[i]+=mu[i-1];
}

ll solve(ll a,ll b,ll c,ll d)
{
    ll ret=0;
    for(ll l=1,r;l<=b;l=r+1)
    {
        r=min(b/(b/l),d/(d/l));
        if((a-1)/l) r=min(r,(a-1)/((a-1)/l));
        if((c-1)/l) r=min(r,(c-1)/((c-1)/l));
        ret+=(mu[r]-mu[l-1])*(b/l-(a-1)/l)*(d/l-(c-1)/l);
    }
    return ret;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0); make_prime();
    int t; cin>>t;
    while(t--)
    {
        ll a,b,c,d,k; cin>>a>>b>>c>>d>>k;
        if(b>d) swap(a,c),swap(b,d);
        cout<<solve(ceil(a,k),b/k,ceil(c,k),d/k)<<endl;
    }
    return 0;
}