洛谷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;
}