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

题意：

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

Solution：

即计算
∑ a 1 = l r ∑ a 2 = l r . . . ∑ a n = l r [ g c d ( a 1 , a 2 , . . . , a n ) = k ] \sum_{a_{1}=l}^{r}\sum_{a_{2}=l}^{r}...\sum_{a_{n}=l}^{r}[gcd(a_{1},a_{2},...,a_{n})=k] a1=l∑ra2=l∑r...an=l∑r[gcd(a1,a2,...,an)=k]
莫比乌斯反演即计算
∑ a 1 = ⌈ l k ⌉ ⌊ r k ⌋ ∑ a 2 = ⌈ l k ⌉ ⌊ r k ⌋ . . . ∑ a n = ⌈ l k ⌉ ⌊ r k ⌋ [ g c d ( a 1 , a 2 , . . . , a n ) = 1 ] = ∑ a 1 = ⌈ l k ⌉ ⌊ r k ⌋ ∑ a 2 = ⌈ l k ⌉ ⌊ r k ⌋ . . . ∑ a n = ⌈ l k ⌉ ⌊ r k ⌋ ∑ d ∣ g c d ( a 1 , a 2 , . . . , a n ) μ ( d ) \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) a1=⌈kl⌉∑⌊kr⌋a2=⌈kl⌉∑⌊kr⌋...an=⌈kl⌉∑⌊kr⌋[gcd(a1,a2,...,an)=1]=a1=⌈kl⌉∑⌊kr⌋a2=⌈kl⌉∑⌊kr⌋...an=⌈kl⌉∑⌊kr⌋d∣gcd(a1,a2,...,an)∑μ(d)
设 L = ⌈ l k ⌉ , R = ⌊ r k ⌋ L=\lceil \frac{l}{k} \rceil,R=\lfloor \frac{r}{k} \rfloor L=⌈kl⌉,R=⌊kr⌋，优先枚举因子，即计算
∑ d = 1 R μ ( d ) ( ⌊ R d ⌋ − ⌊ L − 1 d ⌋ ) n \sum_{d=1}^{R}\mu(d)(\lfloor \frac{R}{d} \rfloor-\lfloor \frac{L-1}{d} \rfloor)^{n} d=1∑Rμ(d)(⌊dR⌋−⌊dL−1⌋)n
整除分块即可完成此式操作，其中 μ ( x ) \mu(x) μ(x)需要区间回答，需要处理成前缀和相减，前缀和只需要用杜教筛计算

时间复杂度：

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

// #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;
}

stdforces