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

题目链接:
(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)=\sum_{n|d} f(d), f(n)=\sum_{n|d}\mu (\frac{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)\le 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)=a^N-a$$ 其中a为[l,r]内是x的倍数的数的个数。公式解释: 如果是随意选，共有 $a^N$种选法，然后去掉全部相等的选法，选N个全部相等的数就相当于只选一个数，因此有a种选法，从 $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;
}