BZOJ 3930 Luogu P3172 选数 (莫比乌斯反演)

Luogu 同时被 3 个专栏收录
9 篇文章 0 订阅
12 篇文章 0 订阅
2 篇文章 0 订阅

题目链接:
(Luogu)https://www.luogu.org/problemnew/show/P3172
(BZOJ)http://www.lydsy.com/JudgeOnline/problem.php?id=3930

题目大意:
给定N,M,L,R,从区间[L,R]内选出N个整数使得它们的gcd恰好为m,求合法的选数方案数对1e9+7取模的值。1<=N,M,L,R<=1e9, R-L<=1e5.

思路分析:
gcd? 那就莫比乌斯反演好了。
令F(m)表示从[L,R]中选出N个数,其gcd为m的倍数的方案数。
f(m)表示从[L,R]中选出N个数,其gcd 恰好为m方案数。(莫比乌斯反演常见做法)
我们要求的是f(m),为了简化运算,我们令l等于大于等于L的最小的m的倍数,r等于小于等于L的最大的m的倍数。然后l/=m,r/=m,问题转化为求f(1). (莫比乌斯反演常见做法)
根据莫比乌斯反演公式

F(n)=n|df(d),f(n)=n|dμ(dn)F(d) F ( n ) = ∑ n | d f ( d ) , f ( n ) = ∑ n | d μ ( d n ) F ( d )
, F(n)可以O(1)求得,直接反演即可。
现在面临两个问题:
1. F(x)和f(x)的定义域是什么?
2. 如何O(1)求F(x)?

先来解决第二个问题:
F(x)其实就是[l,r]内是x的倍数的数的个数的N次方,可以用快速幂求得。具体见代码getF函数。

难点在于第一个问题:
首先我们知道,定义域不超过r. 而r=R/M是1e9级别的,因此必须优化,发现更多的性质。
F(x)既然表示选出N个数gcd为x的方案数,那我们观察以下式子

gcd(x,y)yx(x<y) gcd ( x , y ) ≤ y − x ( x < y )
如果选的数不全相等,那它们的gcd一定不会超过r-l, 也就是F(x)和f(x)的定义域就会缩小到r-l, 而r-l是1e5级别的!这就很美妙了!
现在只要处理一下选出的所有数全相等的情况了。
为了缩小定义域,我们给F(x)和f(x)分别添加一个条件: F(x)表示表示从[L,R]中选出 不全相等的 N个数,其gcd为 x的倍数的方案数,f(x)表示表示从[L,R]中选出 不全相等的 N个数,其gcd 恰好为x的方案数,枚举定义域[1,r-l]莫比乌斯反演求出f(1)即可。
而定义变了以后,O(1)计算F(x)的方法也出现了变动:
F(x)=aNa F ( x ) = a N − a
其中a为[l,r]内是x的倍数的数的个数。公式解释: 如果是随意选,共有 aN a N 种选法,然后去掉全部相等的选法,选N个全部相等的数就相当于只选一个数,因此有a种选法,从 aN a N 中扣除。
以上是计算f(1)的方法。
f(1)算完后,还要加上从[l,r]中选N个全相等的数使得gcd为1的方案数。那显然唯一方案就是全选1,如果1被包含在区间[l,r]中答案就是f(1)+1,否则答案为f(1).

代码实现:

#include<cstdio>
using namespace std;

const int N = 1e5+1;
const long long P = 1e9+7;
long long n,m,lb,rb;
int mu[N+4];
long long p[N+4];
bool f[N+4];
int pn;

void Mobius()
{
    mu[1] = 1; pn = 0;
    for(int i=2; i<=N; i++)
    {
        if(!f[i]) {pn++; p[pn] = i; mu[i] = -1;}
        for(int j=1; j<=pn && i*p[j]<=N; j++)
        {
            f[p[j]*i] = true;
            if(i%p[j]==0) {mu[i*p[j]] = 0; break;}
            else mu[i*p[j]] = -mu[i];
        }
    }
}

long long quickpow(long long a,long long b)
{
    a %= P;
    long long cur = a,ret = 1ll;
    for(int i=0; b; i++)
    {
        if(b&(1ll<<i)) {ret *= cur; ret %= P; b-=(1ll<<i);}
        cur *= cur; cur %= P;
    }
    return ret;
}

long long getF(long long a)
{
    long long lt,rt;
    if(lb%a>0ll) lt = lb/a+1;
    else lt = lb/a;
    rt = rb/a;
    return (quickpow(rt-lt+1,n)-(rt-lt+1)+P)%P;
}

int main()
{
    Mobius();
    scanf("%lld%lld%lld%lld",&n,&m,&lb,&rb);
    if(lb%m>0ll) lb = lb/m+1;
    else lb = lb/m;
    rb/=m;
    long long nn = rb-lb,ans = 0ll;
    for(int i=1; i<=nn; i++)
    {
        ans += mu[i]*getF(i);
        ans = (ans+P)%P;
    }
    if(lb<=1 && 1<=rb) {ans++; ans%=P;}
    printf("%lld\n",ans);
    return 0;
}
  • 1
    点赞
  • 0
    评论
  • 0
    收藏
  • 打赏
    打赏
  • 扫一扫,分享海报

参与评论 您还未登录,请先 登录 后发表或查看评论
©️2022 CSDN 皮肤主题:编程工作室 设计师:CSDN官方博客 返回首页

打赏作者

suncongbo

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值