洛谷3172：选数

洛谷3172：选数

题意：

从区间\([L,R]\)选取\(N\)个整数，总共有\((R-L+1)^N\)中方案。

问所有方案中，方案中所有数字的\(gcd\)等于\(K\)的方案有多少个。

对结果\(mod\ 1e9+7\)。

数据范围\(1\leq N,K\leq10^9,1\leq L\leq R\leq10^9,R-l\leq 10^5\)。

思路：

设\(F(n)\)表示从\([L,R]\)中选\(N\)个数的\(gcd\)为\(n\)或\(n\)的倍数的情况数。

设\(f(n)\)表示从\([L,R]\)中选\(N\)个数的\(gcd\)为\(n\)的情况数。

很明显有：

\[F(n)=\sum_{n|d}f(d) \]

莫比乌斯反演：

\[f(n)=\sum_{n|d}F(d)\mu(\frac{d}{n}) \]

答案是\(f(K)\)。

但是\(f(K)\)不太好求，因为这个倍数其实不是很方便枚举，\(f(1)\)是相对好求的。

所以我们只需要在区间\([\frac{L}{K},\frac{R}{K}]\)找\(f(1)\)就好了。

但这里有个问题，就是如果\(L\%K\neq 0\)，那么就会有\(\lfloor\frac{L}{K}\rfloor *K<L\)，会使一些不在\([L,R]\)范围内的数混入答案当中，所以这里加个特判，如果\(L\%k\neq 0\)，就有\(L=L/K+1\)。

接着考虑\(F(n)\)怎么求？

我们知道

\[gcd(a,b)=d\\ \]

有：

\[d|a,d|b \]

所以\(F(n)\)就比较好求了，其实就是\([L,R]\)区间范围内，\(n\)的倍数的个数的\(N\)次方。

\[F(n)=(\lfloor\frac{R}{n}\rfloor-\lfloor\frac{L-1}{n}\rfloor )^N \]

答案就是：

\[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 \]

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 5e6;
const int mod = 1e9 + 7;

int primes[maxn+10], cnt;
bool vis[maxn+10];
ll mu[maxn+10];
void init(int n)
{
    mu[1] = 1;
    for(int i = 2; i <= n; i++)
    {
        if(!vis[i])
        {
            primes[++cnt] = i;
            mu[i] = -1;
        }
        for(int j = 1; primes[j] <= n/i; j++)
        {
            vis[primes[j]*i] = 1;
            if(i % primes[j] == 0) break;
            else mu[i*primes[j]] = -mu[i];
        }
    }
    for(int i = 1; i <= n; i++)
        mu[i] += mu[i-1];
}

ll qmi(ll a, ll b)
{
    a %= mod;
    ll res = 1; res %= mod;
    while(b)
    {
        if(b&1) res = (res*a)%mod;
        a = (a*a)%mod;
        b >>= 1;
    } return res%mod;
}

unordered_map<ll, ll> Smu;
inline ll getSmu(ll n)
{
    if(n <= maxn) return mu[n];
    if(Smu[n]) return Smu[n];
    ll res = 1;
    for(ll l = 2, r; l <= n; l = r+1)
    {
        r = n/(n/l);
        res -= (r-l+1)%mod*getSmu(n/l)%mod;
        res = (res%mod+mod)%mod;
    } return Smu[n] = res;
}

ll N, K, L, R;

int main()
{
    init(maxn);
    cin >> N >> K >> L >> R;
    R /= K;
    if(L % K == 0) L = L/K;
    else L = L/K+1;

    L = L-1; ll ans = 0;
    for(ll l = 1, r; l <= R; l = r+1)
    {
        if(L/l == 0) r = R/(R/l);
        else r = min(R/(R/l), L/(L/l));
        ans = (ans+qmi(R/l-L/l, N)%mod*(getSmu(r)-getSmu(l-1))%mod)%mod;
    }
    cout << (ans%mod+mod)%mod << endl;
    return 0;
}