跳转至

数论分块（整除分块）

取整函数引理

上下取整函数有一些基本的不等式转化：

\[ \lfloor x \rfloor \ge k \Longleftrightarrow x \ge k\\ \lfloor x \rfloor < k \Longleftrightarrow x < k\\ \lceil x \rceil > k \Longleftrightarrow x > k\\ \lceil x \rceil \le k \Longleftrightarrow x \le k\\ \]

数论分块

数论分块可以在 \(O(\sqrt n)\) 的复杂度内求出：

\[ \sum_{i=1}^{m}{f\left(\lfloor\frac{n}{i}\rfloor\right)} \]

注意到，\(\lfloor \frac{n}{i} \rfloor\) 最多只有 \(2\sqrt{n}\) 个值（根号分治）。我们考虑将相同的值一次计算出来。最重要的问题就是知道区间的起点和终点，即

\[ \begin{align*} l&=\min_{\lfloor \frac{n}{i}\rfloor = k}{\{i\}}\\ r&=\max_{\lfloor \frac{n}{i}\rfloor = k}{\{i\}} \end{align*} \]

我们用 \(l_i\) 和 \(r_i\) 表示第 \(i\) 段区间的左右边界。容易发现，对于 \(i>1\)，\(l_i=r_{i-1}+1\)。又因为 \(l_1=1\)，我们只要知道前一个区间的右边界，就能知道当前的左边界了。

假设我们已知 \(l_i\)，要求出 \(r_i\)。设 \(k=\lfloor \frac{n}{l_i} \rfloor\)，\(r_i\) 为使 \(\lfloor \frac{n}{r_i} \rfloor \ge k\) 的最大值。我们注意到 \(\lfloor \frac{n}{r_i} \rfloor \ge k\) 等价于 \(\frac{n}{r_i} \ge k\)，从而

\[ r_i\le \frac{n}{k} \]

因为 \(r_i\) 是整数，因此 \(r_i = \lfloor \frac{n}{k}\rfloor\)，即

\[ r_i=\left\lfloor{\frac{n}{\lfloor \frac{n}{l_i} \rfloor}}\right\rfloor \]

写一个 for 循环，循环求出 \(l_i\) 对应的 \(r_i\)，更新答案，并使 \(l_{i+1}=r_i+1\)。时间复杂度 \(O(\sqrt{n})\)。

P2261 [CQOI2007] 余数求和

模板题。

P6788 「EZEC-3」四月樱花

给定正整数 \(n\) 和质数 \(p\)，求

\[ \left(\prod_{x=1}^n\prod_{y|x}\frac{y^{d(y)}}{\prod_{z|y}(z+1)^2}\right)\bmod p \]

其中 \(d(y)\) 表示 \(y\) 的因数个数。

\(n\le 2.5\times 10^9\)

先化简式子。注意到 \(z\) 枚举的是 \(y\) 的因数，因此可以把 \(y^{d(y)}\) 拆到第三层 \(\prod\) 中：

\[ \prod_{x=1}^n\prod_{y|x}\frac{y^{d(y)}}{\prod_{z|y}(z+1)^2}=\prod_{x=1}^n\prod_{y|x}\prod_{z|y}\frac{y}{(z+1)^2} \]

后面的式子和 \(y,z\) 同时有关，还是不好看。由于 \(z\) 取遍 \(y\) 的全体因数，因此考虑将 \(y\) 拆成 \(z\cdot y/z\)：

\[ \begin{align*} \prod_{x=1}^n\prod_{y|x}\prod_{z|y}\frac{y}{(z+1)^2}&=\prod_{x=1}^n\prod_{y|x}\prod_{z|y}\frac{z\cdot y/z}{(z+1)^2}\\ &=\prod_{x=1}^n\prod_{y|x}\prod_{z|y}\frac{z^2}{(z+1)^2} \end{align*} \]

现在后面的式子只和 \(z\) 有关了。考虑枚举 \(z\)：

\[ \begin{align*} \prod_{x=1}^n\prod_{y|x}\prod_{z|y}\frac{z^2}{(z+1)^2}&=\prod_{z=1}^n\left(\frac{z^2}{(z+1)^2}\right)^{\small\sum\limits_{z|y}\sum\limits_{y|x}1}\\ &=\prod_{z=1}^n\left(\frac{z^2}{(z+1)^2}\right)^{\small\sum\limits_{y=1}^{n/z}\sum\limits_{x=1}^{n/yz}1}\\ &=\prod_{z=1}^n\left(\frac{z^2}{(z+1)^2}\right)^{\small\sum\limits_{y=1}^{n/z}\lfloor\frac{n}{yz}\rfloor}\\ &=\prod_{z=1}^n\left(\frac{z^2}{(z+1)^2}\right)^{\small\sum\limits_{y=1}^{n/z}\big\lfloor\frac{\lfloor\frac{n}{z}\rfloor}{y}\big\rfloor} \end{align*} \]

设 \(f(x)=\sum_{i=1}^{x}\lfloor\frac xi\rfloor\)，那么上式又可以写成：

\[ \prod_{z=1}^n\left(\frac{z^2}{(z+1)^2}\right)^{\small\sum\limits_{y=1}^{n/z}\big\lfloor\frac{\lfloor\frac{n}{z}\rfloor}{y}\big\rfloor}=\prod_{z=1}^n\left(\frac{z^2}{(z+1)^2}\right)^{\small f(\lfloor\frac{n}{z}\rfloor)} \]

考虑对 \(z\) 整除分块，发现 \(z\) 取一段区间时，底数部分会约分成 \(\frac{l^2}{(r+1)^2}\)，而指数部分我们可以再套一层整除分块。这样，时间复杂度为

\[ O\Biggl(\sum_{z=1}^{\sqrt n}\sqrt z+\sqrt\frac{n}{z}\Bigg)=O(n^{3/4}) \]

其实可以进一步优化。考虑预处理出 \(f\) 的一段前缀，这样可以减少后面的复杂度。考虑经典转化：

\[ \sum_{i=1}^{n}\lfloor\frac{n}{i}\rfloor=\sum_{1\le i\le j\le n}[i|j]=\sum_{i=1}^{n}d(i) \]

也就是因数个数函数的前缀和。由于因数个数函数 \(d(x)\) 是积性函数，因此可以用线性筛在 \(O(B)\) 的时间复杂度内得到 \(f\) 的前 \(B\) 项。考虑平衡时间复杂度，预处理前 \(B\) 项的时间复杂度为

\[ O\Biggl(B+\sum_{i=1}^{n/B}\sqrt\frac ni\Bigg)=O(B+\frac{n}{\sqrt B}) \]

在 \(B=n^{2/3}\approx10^6\) 时最小，为 \(O(n^{2/3})\)。算上求逆元，时间复杂度多带一个 \(\log\)。

代码 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#include<iostream>
#define int long long
using namespace std;
const int N2 = 2e6 + 10;

int n, MOD, ans;
int npri[N2], pri[N2], c[N2], pre[N2], d[N2], pcnt;

int s[N2];

int f(int k) {
    if(k < N2) return s[k];
    int res = 0;
    for(int l = 1, r; l <= k; l = r + 1) {
        r = min(k, k / (k / l));
        res = (res + (r - l + 1) * (k / l) % (MOD - 1)) % (MOD - 1);
    }
    return res;
}

inline int pw2(int x) { return x * x % MOD; }
inline int qpow(int a, int b) {
    int res = 1;
    while(b) {
        if(b & 1) res = res * a % MOD;
        a = a * a % MOD;
        b >>= 1;
    }
    return res;
}
inline int inv(int a) { return qpow(a, MOD - 2); }

signed main() {

    cin >> n >> MOD;

    d[1] = 1; c[1] = 0; 
    for(int i = 2; i < N2; i++) {
        if(!npri[i]) { pri[++pcnt] = i; d[i] = 2; c[i] = 1; pre[i] = 1; }
        for(int j = 1; j <= pcnt; j++) {
            if(i * pri[j] >= N2) break;
            npri[i * pri[j]] = 1;
            if(i % pri[j] == 0) {
                c[i * pri[j]] = c[i] + 1; pre[i * pri[j]] = pre[i];
                d[i * pri[j]] = d[pre[i]] * (c[i] + 2) % (MOD - 1);
                break;
            } else {
                c[i * pri[j]] = 1; pre[i * pri[j]] = i;
                d[i * pri[j]] = d[i] * 2 % (MOD - 1);
            }
        }
    }

    for(int i = 1; i < N2; i++) s[i] = (s[i - 1] + d[i]) % (MOD - 1);

    ans = 1;
    for(int l = 1, r; l <= n; l = r + 1) {
        r = min(n, n / (n / l));
        ans = ans * qpow(pw2(l) * pw2(inv(r + 1)) % MOD, f(n / l)) % MOD;
    }

    cout << ans << '\n';

    return 0;
}