牛客P21546 莫比乌斯反演+杜教筛

最新推荐文章于 2025-11-24 15:27:16 发布

原创最新推荐文章于 2025-11-24 15:27:16 发布 · 299 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#算法

题意：

给出 $n, k, l, r$ ，从区间 $[l, r]$ 内取出 $n$ 个数，并且他们的最大公约数为 $k$ ，有多少种取法？这 $n$ 个数可以有相等的

Solution：

即计算
$\sum_{a_{1}=l}^{r}\sum_{a_{2}=l}^{r}...\sum_{a_{n}=l}^{r}[gcd(a_{1},a_{2},...,a_{n})=k]$
莫比乌斯反演即计算
$\sum_{a_{1}=\lceil \frac{l}{k} \rceil}^{\lfloor \frac{r}{k} \rfloor}\sum_{a_{2}=\lceil \frac{l}{k} \rceil}^{\lfloor \frac{r}{k} \rfloor}...\sum_{a_{n}=\lceil \frac{l}{k} \rceil}^{\lfloor \frac{r}{k} \rfloor}[gcd(a_{1},a_{2},...,a_{n})=1]=\sum_{a_{1}=\lceil \frac{l}{k} \rceil}^{\lfloor \frac{r}{k} \rfloor}\sum_{a_{2}=\lceil \frac{l}{k} \rceil}^{\lfloor \frac{r}{k} \rfloor}...\sum_{a_{n}=\lceil \frac{l}{k} \rceil}^{\lfloor \frac{r}{k} \rfloor}\sum_{d|gcd(a_{1},a_{2},...,a_{n})}\mu(d)$
设 $L=⌈lk⌉,R=⌊rk⌋L=\lceil \frac{l}{k} \rceil,R=\lfloor \frac{r}{k} \rfloor$ ，优先枚举因子，即计算
$\sum_{d=1}^{R}\mu(d)(\lfloor \frac{R}{d} \rfloor-\lfloor \frac{L-1}{d} \rfloor)^{n}$
整除分块即可完成此式操作，其中 $μ(x)\mu(x)$ 需要区间回答，需要处理成前缀和相减，前缀和只需要用杜教筛计算

时间复杂度：

最坏的情况下，整除分块时计算前缀和会递归计算所有 $i∈[1,n]i\in[1,n]$ 的前缀和，每次杜教筛可以计算 $n\sqrt{n}$ 个前缀和，于是最多 $n\sqrt{n}$ 次杜教筛，因此最坏复杂度为 $O(nn23)=O(n76)O(\sqrt{n}n^{\frac{2}{3}})=O(n^{\frac{7}{6}})$ ，实际复杂度远小于这个复杂度

// #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=1e6+5,inf=0x3fffffff;
const long long INF=0x3fffffffffffffff,mod=1e9+7;

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

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

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

ll qpow(ll a,ll b)
{
    ll ret=1,base=a;
    while(b)
    {
        if(b&1) ret=ret*base%mod;
        base=base*base%mod;
        b>>=1;
    }
    return ret;
}

map<int,int>map1;

int calcmu(int n)
{
    if(n<N) return mu[n];
    if(map1.count(n)) return map1[n];
    int ret=1;
    for(int l=2,r;l<=n;l=r+1)
    {
        r=min(n,n/(n/l));
        ret=((1ll*ret-1ll*(r-l+1)*calcmu(n/l)%mod)%mod+mod)%mod;
    }
    return map1[n]=ret;
}

int solve(int sta,int en,int n)
{
    int ret=0;
    for(int l=1,r;l<=en;l=r+1)
    {
        r=min(en,en/(en/l));
        if((sta-1)/l) r=min(r,(sta-1)/((sta-1)/l));
        ret=(1ll*ret+qpow(en/l-(sta-1)/l,n)*(((calcmu(r)-calcmu(l-1))%mod+mod)%mod)%mod)%mod;
    }
    return ret;
}

int main()
{
    make_prime();
    int n,k,l,r; cin>>n>>k>>l>>r;
    cout<<solve(ceil(l,k),r/k,n);
    return 0;
}