题目链接：LOJ 2095

我们知道，从区间 $[L, H]$（$L$ 和 $H$ 为整数）中选取 $N$ 个整数，总共有 $(H - L + 1) ^ N$ 种方案。小 Z 很好奇这样选出的数的最大公约数的规律，他决定对每种方案选出的 $N$ 个整数都求一次最大公约数，以便进一步研究。然而他很快发现工作量太大了，于是向你寻求帮助。你的任务很简单，小 Z 会告诉你一个整数 $K$，你需要回答他最大公约数刚好为 $K$ 的选取方案有多少个。由于方案数较大，你只需要输出其除以 $10 ^ 9 + 7$ 的余数即可。

数据范围：$1 \le N, K \le 10 ^ 9​$，$1 \le L \le H \le 10 ^ 9​$，$H - L \le 10 ^ 5​$。

Solution

首先可以轻松地转化为区间 $[L', H']$ 中选择 $n$ 个数使得 $\gcd$ 为 $1$ 的方案数。

我们设 $f(x)$ 为 $\gcd$ 为 $x$ 的方案数，发现无法直接求 $f(1)$，故考虑设 $g(x)$ 为 $\gcd$ 为 $x$ 的倍数的方案数，很显然有：

$$ g(x) = \left(\left\lfloor\frac{H'}{x}\right\rfloor - \left\lfloor\frac{L' - 1}{x}\right\rfloor\right) ^ n = \sum_{x \mid d} f(d) $$

根据莫比乌斯反演，可以得到：

$$ f(x) = \sum_{x \mid d} g(d) \cdot \mu(\frac{d}{x}) $$

由于我们要求的是 $f(1)$，那么：

$$ f(1) = \sum_{i = 1} ^ {H'} g(i) \cdot \mu(i) $$

对 $g(i)$ 我们可以数论分块，但是需要求出 $\mu$ 的前缀和，因此需要杜教筛来解决。

听说有 $\mathcal O((H - L) \log (H - L))$ 的做法？

时间复杂度：$\mathcal O(H ^ {\frac{2}{3}})$。

Code

#include <cstdio>
#include <unordered_map>

const int N = 1e6 + 5, LIM = 1e6;
const int P = 1e9 + 7;

int k, d, l, r, tot, p[N / 10], mu[N], s[N];
bool flg[N];
std::unordered_map<int, int> S;

void add(int &x, int y) {
    (x += y) >= P && (x -= P);
}
void sub(int &x, int y) {
    (x -= y) < 0 && (x += P);
}
int pow(int x, int p) {
    int ans = 1;
    for (; p; p >>= 1, x = 1LL * x * x % P) {
        if (p & 1) ans = 1LL * ans * x % P;
    }
    return ans;
}
void sieve(int n) {
    std::fill(flg + 1, flg + n + 1, true);
    mu[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (flg[i]) {
            p[++tot] = i, mu[i] = -1;
        }
        for (int j = 1; j <= tot && i * p[j] <= n; j++) {
            flg[i * p[j]] = false;
            if (i % p[j] == 0) {
                mu[i * p[j]] = 0;
                break;
            } else {
                mu[i * p[j]] = -mu[i];
            }
        }
    }
    for (int i = 1; i <= n; i++) {
        if (mu[i] < 0) mu[i] += P;
        add(s[i] = mu[i], s[i - 1]);
    }
}
int F(int n) {
    if (n <= LIM) return s[n];
    if (S.find(n) != S.end()) return S[n];
    int ans = (n >= 1);
    for (int l = 2, r; l <= n; l = r + 1) {
        r = n / (n / l);
        sub(ans, 1LL * (r - l + 1) * F(n / l) % P);
    }
    return S[n] = ans;
}
int solve(int n, int L, int R) {
    int ans = 0;
    for (int l = 1, r; l <= n; l = r + 1) {
        r = std::min(L < l ? n : L / (L / l), R / (R / l));
        add(ans, 1LL * pow(R / l - L / l, k) * (F(r) - F(l - 1) + P) % P);
    }
    return ans;
}
int main() {
    sieve(LIM);
    scanf("%d%d%d%d", &k, &d, &l, &r);
    l = (l - 1) / d + 1, r /= d;
    printf("%d\n", solve(r, l - 1, r));
    return 0;
}