数论算法学习笔记

· · 算法·理论

初级算法（知识点）

组合数

Lucas 定理

对于

n,m\in\mathbb{N}^+,p\in prime

若 n,m 在 p 进制下表示为

n=\overline{a_1a_2\ldots a_k} m=\overline{b_1b_2\ldots b_k}

其中，如果不足 k 位则补 0。则

C_n^m=\prod_{i=1}^kC_{a_i}^{b_i}\mod p

证明

我们先对 n,m 归纳。假设当前对 \forall\ 1\le i\lt n,1\le j\lt m 均成立。

令

n=sp+q,m=tp+r,0\le q\lt p,0\le r\lt p

我们构造一个

\begin{aligned}f(x)&=(1+x)^n\\&=(1+x)^{sp+q}\\&=((1+x)^p)^s(1+x)^q\\&\equiv(1+x^p)^s(1+x)^q\\&=(C_s^0+C_s^1x^p+C_s^2x^{2p}+\ldots+C_s^sx^{sp})(C_q^0+C_q^1x^1+\ldots+C_q^qx^q)\mod p\end{aligned}

由此可见，因为 0\le q\lt p，所以 x_m=x^{tp+r} 的系数必然是 x^{tp} 在 (1+x^p)^s 中的系数与 x^r 在 (1+x)^q 中的系数之积，也就是

C_s^t\times C_q^r=C_n^m

归纳成立。

一些疑惑：

题目

中国剩余定理

对于

a_1,a_2,\ldots,a_n,b_1,b_2,\ldots,b_n\in\mathbb{N}^+\text{ 且满足 }a_1,a_2,\ldots,a_n\text{ 两两互质}

则方程

\left\{\begin{matrix}x\equiv b_1\mod a_1\\x\equiv b_2\mod a_2\\\vdots\\x\equiv b_n\mod a_n\\\end{matrix}\right.

必然有解。

证明

依然归纳。容易发现，对于前 k 个约束条件，设已有解 x，那么

x+a_1a_2\ldots a_k

依然是解。

由于 \{a\} 中的数两两互质，得到

\{x,x+a_1a_2\ldots a_k,x+2a_1a_2\ldots a_k,\ldots,x+(a_{k+1}-1)a_1a_2\ldots a_k\}

中的数形成 \mod a_{k+1} 的完全剩余系，其中必然有被 a_{k+1} 整除的。

归纳成立。

一些疑惑

题目

裴蜀定理（应用：exgcd）

十分重要的定理（算法）。

对于 x,y\in\mathbb{Z}，\exist\ u,v\in\mathbb{Z} 使得 xu+yv=\gcd(x,y)，并且对 \forall\ u,v\in\mathbb{Z} 必有 \gcd(x,y)\mid(xu+yv)。

证明

我们先记 d=\gcd(x,y)。

首先，第二问是非常容易的，因为 d\mid xu 且 d\mid yv，所以必然有 d\mid(xu+yv)。

对于第一问，设 x=x_1d，y=y_1d，这样有 (x_1,y_1)=1，我们只需证明 \exist\ u,v\in\mathbb{Z} 使得 x_1u+y_1v=1

与中国剩余定理类似，因为 (x_1,y_1)=1，所以 \{0,y_1,2y_1,\ldots,(x_1-1)y_1\} 构成 \mod x_1 的完全剩余系，其中必 \exist\ 0\le v\lt x_1 使得 x_1\mid(1-y_1v)，再令 u=\frac{1-y_1v}{x_1} 即可。

exgcd（拓展欧几里得算法）就是用来求一组满足要求的 (u,v) 的算法。代码模板如下：

int exgcd(int x, int y, int &u, int &v) {
    if (!y) return u = 1, v = 0, x;
    ans = exgcd(y, x % y, u, v);
    int t = u; u = v, v = t - x / y * v;
    return ans;
}

一些疑惑

这次没有疑惑。o(*≧▽≦)ツ

题目

其它内容

如果想要了解某些定理的拓展版本，请前往模板题或 oi-wiki。

高级一点的算法

数论函数

定义：若 f(x) 满足其定义域、值域均为正整数集则称其为数论函数。

进一步地，若数论函数 f(x) 满足对 x,y\in\mathbb{N}^+ 且 (x,y)=1，均满足 f(xy)=f(x)f(y)，则称 f(x) 为积性函数。

常见的数论函数包括 \left\{\begin{matrix}\epsilon(x)=[x=1]\\1(x)=1\\\text{id}(x)=x\\\varphi(x)\\\mu(x)\end{matrix}\right. 等。

其中 \epsilon(x) 为迪利克雷卷积意义下的单位元，\varphi(x)（欧拉函数）表示 \lt x 的数中与 x 互质的数的数量，\mu(x)（莫比乌斯函数）表示 1(x) 关于单位元 \epsilon(x) 的逆。

具体地，若记 x=\prod_{i=1}^kp_i^{\alpha_i}，其中 p_i 均为质数，\alpha_i 均为正整数，则

\varphi(x)=x\times\prod_{i=1}^k\frac{p_i-1}{p_i} \mu(x)=\left\{\begin{matrix}-1,&[2\nmid k]\land[\alpha_i=1]\\0,&\exist\ \alpha_i\ge2\\1,&[2\mid k]\land[\alpha_i=1]\end{matrix}\right.

这里定义两个数论函数 f(x)，g(x) 的迪利克雷卷积为

h(x)=\sum_{k\mid x}f(k)g(\frac xk)

并且称 f 为 g 关于 h 的逆函数。

线性筛

也称欧拉筛，一种可以在 \mathcal{O}(n) 的时间复杂度内求出 n 以内的所有质数的筛法。它也可以做到同时筛出 n 以内所有 \varphi(x) 以及 \mu(x) 的值。其思想就是根据每个数的最小质因子来更新这个数，这样能够确保时间复杂度是线性的。

代码模板如下：

void sieve(int N) {
    phi[1] = miu[1] = 1;
    for (int i = 2; i <= N; i++) {
        if (!is_p[i]) p[++c] = i, phi[i] = i - 1, miu[i] = -1;
        for (int j = 1; j <= c && i * p[j] <= N; j++) {
            is_p[i * p[j]] = 1;
            if (i % p[j] == 0) { // 保证 p[j] 是 i 的最小质因子
                phi[i * p[j]] = phi[i] * p[j];
                miu[i * p[j]] = 0;
                break;
            }
            phi[i * p[j]] = phi[i] * (p[j] - 1);
            miu[i * p[j]] = -miu[i];
        }
    }
}

莫比乌斯反演

因为莫比乌斯函数 \mu(x) 有着 1*\mu=\epsilon 的关键性质，我们可以得到这个：

感性理解一下：我们如果认为 \mu=\frac\epsilon1（当然，这样是不对的，但是理解一下就行了），并且 1=\frac gf，那么 \frac fg=\frac\epsilon1=\mu，故 f=g*\mu。

转化成代数式就是：

\because g(x)=\sum_{k\mid x}f(k) \therefore f(x)=\sum_{k\mid x}\mu(\frac xk)g(k)

但是，更常用的却是另一个性质：

\because g(x)=\sum_{x\mid k}f(k) \therefore f(x)=\sum_{x\mid k}\mu(\frac kx)g(k)

但是，最常用的却是 x=1 的情况：

f(1)=\sum_{k}\mu(k)g(k)

介于 \mu(x) 可以通过各种筛法快速得出，而 g(x) 必然是很好求的函数（不然你为什么会用它呢？），我们就可以快速得出 f(x) 的值。

题目

这些都是很好的练习题（当然，求答案时肯定要用到很多科技）。

整除分块

几乎所有筛法必要的技能。

例如，如果我们要求 \sum_{i=1}^n\lfloor\frac ni\rfloor，正常枚举是 \mathcal O(n) 的，但是我们要一种更快的做法。

为了求这个，我们必须证明，\lfloor\frac ni\rfloor 只有 \mathcal O(\sqrt n) 种取值。

证明

证明完毕。

那么，我们就可以利用这个整除分块了，假设我们可以求得 l,r 使得 \lfloor\frac nl\rfloor=\lfloor\frac nr\rfloor，那么区间 [l,r] 中的取值之和就是 (r-l+1)\times\lfloor\frac nl\rfloor。

具体求 l,r 代码如下：

for (int l = 1, r; l <= n; l = r + 1)
    r = n / (n / l),
    ans += (r - l + 1) * (n / r);

时间复杂度为 \mathcal O(\sqrt n)。

杜教筛

一种比线性筛快那么亿点点，有局限性，但是十分好写好吃的筛法。

为了方便，我们记所有 \lfloor\frac ni\rfloor 为 n 的整除点。

假设我们要求积性函数 g(x) 在 n 所有整除点上前缀和。

首先，为了满足筛法要求我们设 f*g=h，那么，由卷积定义式，我们有

\begin{aligned}\sum_{i=1}^nh(i)&=\sum_{ij\le n}f(i)g(j)\\&=\sum_{i=1}^n\sum_{j=1}^{\lfloor\frac ni\rfloor}f(i)g(j)\\&=\sum_{i=1}^nf(i)G(\lfloor\frac ni\rfloor)\end{aligned}

这里 G(x) 指 g(x) 的前缀和。求的东西就是 G(x) 在整除点上的点值。

假设我们已经求出了 n 的所有整除点（这次不包括自己）上 G 的点值，那么我们只需要求 G(n) 的值。

G(n)=\frac{\sum_{i=1}^nh(i)-\sum_{i=\large\color{red}2}^nf(i)G(\lfloor\frac ni\rfloor)}{f(1)}

而因为 f 是积性函数，必有 f(1)=1，所以我们得到了 G(n)：

G(n)=\sum_{i=1}^nh(i)-\sum_{i=2}^nf(i)G(\lfloor\frac ni\rfloor)

如果 f,h 的前缀和可以快速求出，那么后面的那一坨就可以通过整除分块 \mathcal O(\sqrt n) 求出。

所以，我们只需要从小到大遍历所有整除点，然后对每个整除点都做一遍就可以了！

使用前提

\left\{\begin{matrix}1*\varphi=\text{id}\\1*\mu=\epsilon\end{matrix}\right.

时间复杂度

杜教筛的时间复杂度 \mathcal O(n^{3/4})，但是证明较为复杂。

复杂度瓶颈在于对每个整除点整除分块，依然使用我们证明整除分块时间复杂度为 \mathcal O(\sqrt n) 时的方法，分 i 与 \sqrt n 的关系讨论。

复杂度为

\sum_{i=1}^{\sqrt n}\sqrt{\frac ni}+\sum_{i=1}^{\sqrt n}\sqrt i

显然后面的式子可以直接忽略，所以

\begin{aligned}&\mathcal O(\sum_{i=1}^{\sqrt n}\sqrt{\frac ni}+\sum_{i=1}^{\sqrt n}\sqrt i)\\=&\mathcal O(\sum_{i=1}^{\sqrt n}\sqrt{\frac ni})\\=&\mathcal O(\sqrt n\sum_{i=1}^{\sqrt n}\sqrt{\frac1i})\\=&\mathcal O(\sqrt n\sum_{i=1}^{\sqrt n}i^{-\frac12})\\=&\mathcal O(\sqrt n\int_0^{\sqrt n}x^{-\frac12}dx)\\=&\mathcal O(\sqrt n(\left. x^{\frac12}\right|_0^{\sqrt n}))\\=&\mathcal O(\sqrt n\sqrt{\sqrt n})\\=&\mathcal O(n^{\frac34})\end{aligned}

依然不够优秀。怎么办？想到可以预处理前 m 个 g(x)，这样我们再算一下时间复杂度：

\begin{aligned}&\mathcal O(m+\sqrt n\sum_{i=1}^{\frac nm}i^{-\frac12})\\=&\mathcal O(m+\sqrt n\sqrt{\frac nm})\end{aligned}

（部分过程同上，已省略）

那么，为了平衡复杂度，解方程

m=\sqrt n\sqrt{\frac nm}

得 m=n^{2/3}，故杜教筛的最终时间复杂度是 \mathcal O(n^{2/3})。

当然，介于存储用的 map 以及整除分块大常数，杜教筛没有预想的那么快，不过还是足够了。

题目

有一些题目在莫比乌斯反演中出现过，但是它们也需要用到杜教筛。

例题选讲

题目大意：在区间 [l,r] 中选 n 个数，使得这 n 个数的 \gcd 是 k，问共有多少种取法（一个数可以取多次）。

首先，注意到取的数必然是 k 的倍数，那么其余的数就失去了意义，我们可以直接令 l\leftarrow \lfloor\frac{l-1}k\rfloor，r\leftarrow\lfloor\frac rk\rfloor，这样问题转化为在 (l,r] 中选 n 个数使得它们互质。

我们记 f(x) 为在 (l,r] 中取 n 个数使得它们的最大公约数是 x 的方案数，g(x) 为在 (l,r] 中取 n 个数使得它们的最大公约数是 x 的倍数的方案数。那么

\left\{\begin{matrix}g(x)=(\lfloor\frac rx\rfloor-\lfloor\frac lx\rfloor)^n\\g(x)=\sum_{x\mid k}f(k)\end{matrix}\right.

于是使用莫比乌斯反演，

f(x)=\sum_{x\mid k}\mu(\frac kx)g(k)

在 x=1 时有

\begin{aligned}f(1)&=\sum_{k}\mu(k)g(k)\\&=\sum_{k}\mu(k)(\lfloor\frac rk\rfloor-\lfloor\frac lk\rfloor)^n\end{aligned}

杜教筛+整除分块处理即可。

#include <bits/stdc++.h>
#define int long long
using namespace std;
const int mod = 1e9 + 7;
const int N = 5e6 + 5;
int miu[N], p[N];
bool is_p[N];
int n, k, L, R, c, ans;
map <int, int> mp;
void sieve(int N) { // 线性筛
    miu[1] = 1;
    for (int i = 2; i <= N; i++) {
        if (!is_p[i]) p[++c] = i, miu[i] = -1;
        for (int j = 1; j <= c && i * p[j] <= N; j++) {
            is_p[i * p[j]] = 1;
            if (i % p[j] == 0) {
                miu[i * p[j]] = 0;
                break;
            }
            miu[i * p[j]] = -miu[i];
        }
    }
    for (int i = 2; i <= N; i++)
        miu[i] += miu[i - 1];
}
int qpow(int x, int y, int res = 1) {
    while (y) {
        if (y & 1) (res *= x) %= mod;
        (x *= x) %= mod; y >>= 1;
    }
    return res;
}
int cal_miu(int x) { // 杜教筛求 miu[x] 前缀和
    if (x <= 5e6) return miu[x];
    if (mp.count(x)) return mp[x];
    int res = 1;
    for (int l = 2, r; l <= x; l = r + 1)
        r = x / (x / l),
        (res -= (r - l + 1) * cal_miu(x / r)) %= mod;
    return mp[x] = res;
}
signed main() {
    scanf("%lld%lld%lld%lld", &n, &k, &L, &R);
    L = (L - 1) / k, R /= k; sieve(5e6);
    for (int l = 1, r; l <= R; l = r + 1) { // 整除分块（注意除以 0！）
        r = R / (R / l); if (L / l) r = min(r, L / (L / l));
        (ans += (cal_miu(r) - cal_miu(l - 1)) * qpow(R / r - L / r, n)) %= mod;
    }
    return printf("%lld", (ans + mod) % mod), 0;
}

Powerful Number 筛（PN 筛）

有时我们可能遇到杜教筛无法解决的函数，这时我们就可能需要用到这个筛法。

首先，我们定义 Powerful Number 为形如 a^2b^3 的数，或者，如果一个数满足它的所有质因子均至少二次，那么我们就称其为 Powerful Number。

我们现在证明 n 以内的 Powerful Number 是 \mathcal O(\sqrt n) 级别的。我们只需枚举 a：

\begin{aligned}&\mathcal O(\sum_{a=1}^{\sqrt n}\sqrt[3]{\frac n{a^2}})\\=&\mathcal O(\sqrt[3]n\sum_{a=1}^{\sqrt n}a^{-\frac23})\\=&\mathcal O(\sqrt[3]n\int_0^{\sqrt n}x^{-\frac23}dx)\\=&\mathcal O(\sqrt[3]n\sqrt{\sqrt[3]n})\\=&\mathcal O(\sqrt n)\end{aligned}

证明完毕。

假设我们要求的前缀和函数为 f(x)。那么我们先找到一个 g(x) 满足 g 在所有质数上的点值和 f 相同。此时，我们再构造一个 h(x) 使得 f=g*h。

首先，由于 h 是积性函数，h(1)=1。并且，对于 \forall\ p 为质数，我们有 f(p)=g(p)h(1)+g(1)h(p)，又因为 f(p)=g(p)，故 h(p)=0。

再次利用 h 是积性函数，如果 x 不是 Powerful Number，那么 x 必有一个质因子次数为 1，假设是 p，那么 h(x)=h(p)h(\frac xp)=0。故 h(x) 只有在 Powerful Number 上的点值可能非零。

这样我们就能算 f 的第 n 项前缀和。

\sum_{i=1}^nf(i)=\sum_{i=1}^nh(i)\sum_{j=1}^{\lfloor\frac ni\rfloor}g(j)

如果我们能够用科技求出 g 在 n 的整除点上的前缀和，那么我们只需要找到 n 以内所有的 Powerful Number，就可以用求 g 的科技的复杂度为瓶颈求出 \sum_{i=1}^nf(i)。

题目

例题选讲

虽然写的是 Min_25 筛，但是我们其实可以用 PN 筛做。

根据筛法条文，我们应该先构造 g(p)=f(p)，p 为质数且 g 在 n 的整除点上的前缀和较容易计算。

注意到 f(p)=p(p-1)=p\varphi(p)，所以我们可以令 g(x)=x\varphi(x)。

接下来我们就需要知道，怎样快速求出 g 整除点上的的前缀和。它虽然并不是简单的 \varphi(x)，但也只是加了一个 x 而已，我们根据古老的公式 1*\varphi=\text{id}，考虑在两边同时乘以 \text{id}，得到 \text{id}*g=\text{id}^2。于是，我们就可以通过杜教筛筛出 g。

最后一个难点就是，如何求出 h=\frac fg 在 Powerful Number 处的点值？

考虑将 Powerful Number 肢解，让它没那么有 Power分解质因数，根据 h 积性函数的特点，我们只需求出 h(p^k)。

并且，因为 f=g*h

\begin{aligned}f(p^k)&=\sum_{x\mid p^k}g(x)h(\frac{p^k}x)\\&=\sum_{i=0}^kg(p^i)h(p^{k-i})\end{aligned}

类似杜教筛，我们再次将 h(p^k) 单独拎出来：

h(p^k)=f(p^k)-\sum_{i=1}^kg(p^i)h(p^{k-i})

但是这样依然不够简洁。我们考虑，能不能直接推出 h(p^k) 的表达式？

先找找规律：

数多了之后渐渐地发现了规律：h(p^k)=(k-1)(p^{k+1}-p^k)，由于太复杂了我不会证qwq

然后根据积性函数的性质求出 h 即可。

求 Powerful Number 的过程用 dfs 就行，枚举当前质数的指数即可。

#include <bits/stdc++.h>
#define int long long
using namespace std;
const int mod = 1e9 + 7;
const int inv2 = 500000004;
const int inv6 = 166666668;
const int N = 5e6 + 5;
int p[1000000], phi[N], c;
bool is_p[N];
void sieve(int N) { // 线性筛
    phi[1] = 1;
    for (int i = 2; i <= N; i++) {
        if (!is_p[i]) p[++c] = i, phi[i] = i - 1;
        for (int j = 1; j <= c && i * p[j] <= N; j++) {
            is_p[i * p[j]] = 1;
            if (i % p[j] == 0) {
                phi[i * p[j]] = phi[i] * p[j];
                break;
            }
            phi[i * p[j]] = phi[i] * (p[j] - 1);
        }
    }
    for (int i = 2; i <= N; i++) (phi[i] *= i) %= mod;
    for (int i = 2; i <= N; i++) (phi[i] += phi[i - 1]) %= mod;
}
int mp[N], n, ans;
int cal_G(int x) { // 杜教筛筛 g
    if (x <= 5e6) return phi[x];
    if (mp[n / x]) return mp[n / x];
    int X = x % mod;
    int res = X * (X + 1) % mod * (2 * X + 1) % mod * inv6 % mod;
    for (int l = 2, r; l <= x; l = r + 1)
        r = x / (x / l),
        (res -= (l + r) % mod * (r - l + 1) % mod * inv2 % mod * cal_G(x / r)) %= mod;
    return mp[n / x] = res;
}
void dfs(int nw, int num, int H) { // num 表示当前 h 值
    if (nw > c || num * p[nw] > n) {
        (ans += H * cal_G(n / num)) %= mod;
        return;
    }
    dfs(nw + 1, num, H); // 当前质数不取
    for (int mu = p[nw] * p[nw], s = 1; num * mu <= n; mu *= p[nw], s++)
        dfs(nw + 1, num * mu, s * (p[nw] - 1) % mod * mu % mod * H % mod); // 计算 h 公式
}
signed main() {
    scanf("%lld", &n), sieve(5e6), dfs(1, 1, 1);
    return printf("%lld", ans), 0;
}

再见~