[CQOI2015] 选数 (莫比乌斯反演 杜教筛)

题意

求从区间 $[L,R]$ 选出 $n$ 个数使得最大公约数为 $k$ 的方案数，对 $10^9+7$ 取模

$1\le n,k\le 10^9$
$1\le L\le R\le 10^9$

分析：

根据题意

$$\sum _{a_1=L}^R\sum _{a_2=L}^R\cdots \sum _{a_n=L}^R[\gcd (a_1,a_2,\cdots ,a_n)=k]$$

把 $k$ 拿到上下界

$$\sum _{a_1=\lfloor \frac{L-1}{k}\rfloor +1}^{\lfloor \frac{R}{k}\rfloor }\sum _{a_2=\lfloor \frac{L-1}{k}\rfloor +1}^{\lfloor \frac{R}{k}\rfloor }\cdots \sum _{a_n=\lfloor \frac{L-1}{k}\rfloor +1}^{\lfloor \frac{R}{k}\rfloor }[\gcd (a_1,a_2,\cdots ,a_n)=1]$$

莫比乌斯反演

$$\sum _{a_1=\lfloor \frac{L-1}{k}\rfloor +1}^{\lfloor \frac{R}{k}\rfloor }\sum _{a_2=\lfloor \frac{L-1}{k}\rfloor +1}^{\lfloor \frac{R}{k}\rfloor }\cdots \sum _{a_n=\lfloor \frac{L-1}{k}\rfloor +1}^{\lfloor \frac{R}{k}\rfloor }\sum _{d\mid \gcd (a_1,a_2,\cdots ,a_n)}\mu (d)$$

交换求和次序

$$\sum _{d=1}^{\lfloor \frac{R}{k}\rfloor }\mu (d)(\lfloor \frac{R}{kd}\rfloor -\lfloor \frac{L-1}{kd}\rfloor {)}^n$$

然后用杜教筛做

代码：

#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N = 2e6 + 5, mod = 1e9 + 7;
int n, k, L, R, primes[N], mobius[N], cnt, res, sum[N];
bool st[N];
unordered_map<int, int> mp;
int qmi(int a, int b) {
    int res = 1;
    while (b) {
        if (b & 1) res = res * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return res;
}
void get_mobius(int n) {
    mobius[1] = 1;
    for (int i = 2; i <= n; i ++) {
        if (!st[i]) {
            primes[cnt ++] = i;
            mobius[i] = -1;
        }
        for (int j = 0; primes[j] * i <= n; j ++) {
            int t = primes[j] * i;
            st[t] = 1;
            if (i % primes[j] == 0) {
                mobius[t] = 0;
                break;
            }
            mobius[t] = -mobius[i];
        }
    }
    for (int i = 1; i <= n; i ++) sum[i] = (sum[i - 1] + mobius[i] + mod) % mod;
}
int Sum(int n) {
    if (n < N) return sum[n];
    if (mp[n]) return mp[n];
    int res = 1;
    for (int l = 2, r; l <= n; l = r + 1) {
        r = n / (n / l);
        res = (res - Sum(n / l) * (r - l + 1) % mod + mod) % mod;
    }
    return mp[n] = res;
}
signed main() {
    get_mobius(N - 1);
    cin >> n >> k >> L >> R;
    L = (L - 1) / k, R /= k;
    for (int l = 1, r; l <= R; l = r + 1) {
        r = min(R / (R / l), L / l ? L / (L / l) : mod);
        res = (res + (Sum(r) - Sum(l - 1)) * qmi(R / l - L / l, n) % mod + mod) % mod;
    }
    cout << res << endl;
}