洛谷P3172 [CQOI2015] 选数题解

题目链接：P3172 [CQOI2015] 选数

题意：

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

由于方案数较大，你只需要输出其除以 \(10^9+7\) 的余数即可。

输入格式：

输入一行，包含四个空格分开的正整数，依次为 \(N,K,L,H\)。

输出格式：

输出一个整数，为所求方案数除以 \(10^9 + 7\) 的余数。

数据范围：

对于 \(100\%\) 的数据，\(1\le N,K\le 10^9\)，\(1\le L\le H\le 10^9\)，\(H-L\le 10^5\)。

设 \(F(k)\) 表示从 \([L,R]\) 中选 \(N\) 个数的最大公约数为 \(k\) 或 \(k\) 的倍数的方案数。

设 \(f(k)\) 表示从 \([L,R]\) 中选 \(N\) 个数的最大公约数是 \(k\) 的方案数。那么 \[ F(k) = \sum_{n \ge k \, \land \, k \,\mid\, n} f(n) \] 莫反，启动！ \[ f(k) =\sum_{n \ge k \, \land \, k \,\mid\, n} \mu(n)F\left(\frac{n}{k}\right) \] 显然有 \[ F(k) = \left(\left\lfloor\frac{R}{k}\right\rfloor-\left\lfloor\frac{L-1}{k}\right\rfloor\right)^N \] 然而即便如此， \(f(K)\) 依然不是很好求，那么我们考虑将 \(L,R\) 均缩小 \(K\) 倍，然后算 \(f(1)\)

这里有一个细节，当 \(K \not\mid L\) 时，\(K\cdot \left\lfloor \frac{L}{K}\right\rfloor < L\) ，可能会多算，所以 \(L\) 要上取整。那么 \[ f(1) = \sum_{i=1}^{\left\lfloor\frac{R}{K}\right\rfloor}\mu(i)F(i) \] 注意上界里是 \(K\) 。那么现在问题就变成了怎么计算 \(\mu\) 的前缀和。

考虑杜教筛，构造 \(g = \mathbf{1}\) ，那么 \((f*g)(n)=\epsilon(n)=[n=1]\) ，都非常好算。

时间复杂度 \(\mathcal{O}(n^{\frac{2}{3}})\) ，其中 \(n\) 表示 \(L,R\) 的值域，并且不受 \(R-L\) 的限制。

代码：

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
const int mod = 1e9 + 7;
void up(int &x, int y) { x < y ? x = y : 0; }
void down(int &x, int y) { x > y ? x = y : 0; }
#define rep(i, a, b) for(int i = (a), i##END = (b); i <= i##END; i++)
#define Rep(i, a, b) for(int i = (a), i##END = (b); i >= i##END; i--)
#define N ((int)(5e6 + 15))

char ck[N]; int pcnt, M, p[N], mu[N], sum[N];
unordered_map<int, int> vis;
void add(int &x, int y) { (x += y) >= mod ? x -= mod : 0; }
int qpow(int a, int b)
{
    int r = 1;
    for(a %= mod; b > 0; b >>= 1, a = a * a % mod)
        if(b & 1) r = r * a % mod;
    return r;
}
void init(const int _ = N - 5)
{
    mu[1] = 1;
    rep(i, 2, _)
    {
        if(!ck[i]) { p[++pcnt] = i; mu[i] = -1; }
        for(int j = 1; j <= pcnt && i * p[j] <= _; j++)
        {
            const int pos = i * p[j]; ck[pos] = true;
            if(i % p[j]) { mu[pos] = -mu[i]; } else { mu[pos] = 0; break; }
        }
    }
    rep(i, 1, _) sum[i] = (sum[i - 1] + mod + mu[i]) % mod;
}
int getsum(int n)
{
    if(n <= M) return sum[n];
    if(vis[n]) return vis[n];
    int res = 1;
    for(int l = 2, r; l <= n; l = r + 1)
    {
        r = n / (n / l);
        add(res, mod - (r + mod - l + 1) % mod * getsum(n / l) % mod);
    }
    return vis[n] = res;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int n, K, L, R; cin >> n >> K >> L >> R; R /= K; L = (L + K - 1) / K;
    init(M = min((int)powl(R, 0.67) + 15, N - 5)); int res = 0; --L;
    for(int l = 1, r; l <= R; l = r + 1)
    {
        if(!(L / l)) r = R / (R / l); else r = min(R / (R / l), L / (L / l));
        add(res, (getsum(r) + mod - getsum(l - 1)) % mod * qpow(R / l - L / l, n) % mod);
    }
    cout << res << '\n';
    return 0;
}