[BZOJ 3930] [CQOI 2015]选数(莫比乌斯反演+杜教筛)

[BZOJ 3930] [CQOI 2015]选数(莫比乌斯反演+杜教筛)

题面

我们知道,从区间\([L,R]\)(L和R为整数)中选取N个整数,总共有\((R-L+1)^N\)种方案。求最大公约数刚好为K的选取方案有多少个。由于方案数较大,你只需要输出其除以1000000007的余数即可。

\[N,K,L,H \leq 10^9,H-L \leq 10^5\]

分析

\(\because \gcd(ka,kb)=k\gcd(a,b)\),我们先把\(L,R\)除以\(K\),然后问题就变成了求gcd=1的方案数

\(f(x)\)表示区间[l,r]里选n个数,gcd为x的方案数

\(F(x)\)表示区间[l,r]里选n个数,gcd被x整除的方案数

\(\because x|\gcd(i,j),\therefore x|i,x|j\)

[l,r]里被x整除的数有\((\lfloor \frac{r}{x} \rfloor-\lfloor \frac{l-1}{x} \rfloor)\)

因此\(F(x)=(\lfloor \frac{r}{x} \rfloor-\lfloor \frac{l-1}{x} \rfloor)^n\)

\(F,f\)显然满足莫比乌斯反演的第二种形式,\(F(x)=\sum_{x|d} f(d)\)

\(f(x)=\sum_{x|d} F(d) \mu(\frac{d}{x})\)

我们要求的是

\[f(1)=\sum_{1|d} F(d) \mu(d)=\sum_{d=1}^r \mu(d) (\lfloor \frac{r}{d} \rfloor-\lfloor \frac{l-1}{d} \rfloor)^n\]

后面的部分可以数论分块然后快速幂求解,但由于\(r \leq 10^9\),不能直接线性筛\(\mu\)的前缀和,需要用杜教筛。


杜教筛:

套路公式:

我们要求\(f\)的前缀和,构造两个函数\(g,h\)满足\(h=f*g\), \(F,G,H\)为它们的前缀和

\[g(1)F(n)=H(n)-\sum_{d=2}^n g(d) F(\frac{n}{d})\]

如果\(f=\mu\),注意到\(\mu*I=\varepsilon\),那么\(g(n)=I(n)=1,h(n)=\varepsilon(n),H(n)=\varepsilon(1)=1\)

代入得\(F(n)=1-\sum_{d=2}^n F(\frac{n}{d})\)

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<map>
#define maxn 2000000
#define mod 1000000007
using namespace std;
typedef long long ll;
int n,k;
ll A,B;

int cnt;
bool vis[maxn+5];
int prime[maxn+5];
int mu[maxn+5];
ll s_mu[maxn+5];
void sieve(int n){
    mu[1]=1;
    for(int i=2;i<=n;i++){
        if(!vis[i]){
            prime[++cnt]=i;
            mu[i]=-1;
        } 
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
            vis[i*prime[j]]=1;
            if(i%prime[j]==0){
                mu[i*prime[j]]=0;
                break;
            }else mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=n;i++) s_mu[i]=(s_mu[i-1]+mu[i])%mod;
}

map<ll,ll>sum_mu;
ll dujiao_sieve(ll x){
    if(x<=maxn) return s_mu[x];
    if(sum_mu.count(x)) return sum_mu[x];
    ll ans=1;
    for(ll l=2,r;l<=x;l=r+1){
        r=x/(x/l);
        ans-=(r-l+1)*dujiao_sieve(x/l)%mod;
        ans=(ans+mod)%mod;
    }
    sum_mu[x]=ans;
    return ans;
} 

inline ll fast_pow(ll x,ll k){
    ll ans=1;
    while(k){
        if(k&1) ans=ans*x%mod;
        x=x*x%mod;
        k>>=1;
    }
    return ans;
}

int main(){
    sieve(maxn);
    scanf("%d %d %lld %lld",&n,&k,&A,&B);
    A=(A-1)/k;
    B/=k;
    ll ans=0;
    for(ll l=1,r;l<=B;l=r+1){
        if(A/l) r=min(A/(A/l),B/(B/l));
        else r=B/(B/l);
//      printf("%d %d\n",l,r);
        ans+=fast_pow(B/l-A/l,n)*(dujiao_sieve(r)-dujiao_sieve(l-1)+mod)%mod;
        ans%=mod;
    }
    printf("%lld\n",ans);
}

转载于:https://www.cnblogs.com/birchtree/p/11438266.html

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值