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

题意：

给出$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$
∑ d = 1 m i n { ⌊ n k ⌋ , ⌊ m k ⌋ } μ ( d ) ∑ i = 1 ⌊ n d k ⌋ ∑ i = 1 ⌊ m d k ⌋ 1 \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 d=1∑min{⌊kn​⌋,⌊km​⌋}​μ(d)i=1∑⌊dkn​⌋​i=1∑⌊dkm​⌋​1
由于$n,m$轮换，不妨设$n\le m$，即
$$\sum _{d=1}^{\lfloor \frac{n}{k}\rfloor }\mu (d)\lfloor \frac{n}{dk}\rfloor \lfloor \frac{m}{dk}\rfloor$$
整除分块即可快速求得此式，$\mu (d)$部分只需要线性筛后处理前缀和，需要注意的是，此处需要求单个$\mu (i)$而不是$\sum \mu (i)$，所以$i\%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\le 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;
}