因此 .org 必須 地址 常見 art names inline esp

手動博客搬家：本文發表於20180310 11:46:11, 原地址https://blog.csdn.net/suncongbo/article/details/79506484

題目鏈接:
(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;
}

BZOJ 3930 Luogu P3172 選數 (莫比烏斯反演)