快速幂与逆元

前言

上一篇 同余与模运算 里，我们学会了在模意义下做加法、减法、乘法。但是乘方呢？比如让你算 $2^{1000000007} \bmod 10^9+7$，你总不能写一个循环乘 $10^9$ 次吧——那得跑 10 秒。

而且上一篇还留了一个坑：模意义下不能直接做除法。$\frac{6}{2} \bmod 5$ 等于 $3$，这没问题。但 $6 \div 3 \bmod 5$ 呢？$\frac{6}{3}=2$，好像也行。那 $\frac{1}{3} \bmod 5$ 呢？小数怎么办？

这篇文章解决这两个问题。而你会发现，它们的答案居然是同一个东西。

问题的本质

快速幂：折半

先忘掉代码，用直觉想一个问题：怎么尽快算出 $2^{16}$？

笨办法：$2 \times 2 \times 2 \times \cdots$，乘 15 次。

聪明办法：

• $2^1 = 2$
• $2^2 = 2 \times 2 = 4$
• $2^4 = 4 \times 4 = 16$
• $2^8 = 16 \times 16 = 256$
• $2^{16} = 256 \times 256 = 65536$

5 次乘法搞定。每次把指数翻倍，$2^{10^{18}}$ 也只需要大约 60 次。这就是”快速幂”——每次把问题规模砍一半。

逆元：模意义下的”倒数”

在普通算术里，$3$ 的倒数是 $\frac{1}{3}$，因为 $3 \times \frac{1}{3} = 1$。

在模 $5$ 的世界里，“倒数”是什么意思？我们要找一个数 $x$，使得 $3 \times x \equiv 1 \pmod{5}$。试一下：$3 \times 1=3$，$3 \times 2=6\equiv1$。所以 $3$ 在模 $5$ 下的”倒数”是 $2$。

这个”倒数”有个正式名字：乘法逆元。费马小定理告诉我们，当模数 $p$ 是素数时，逆元就是 $a^{p-2} \bmod p$——算一次快速幂就得到了。

理论 + 代码

快速幂

把指数 $n$ 写成二进制。比如 $n=13$，二进制是 $1101$，也就是 $n=8+4+1$。所以：

$a^{13} = a^8 \cdot a^4 \cdot a^1$

我们只需要不断自乘得到 $a^1, a^2, a^4, a^8, \ldots$，然后看 $n$ 的哪些二进制位是 1，就把对应的值乘进去。

long long qpow(long long a, long long n, long long mod) {
    long long res = 1;     // res 用来累积最终结果，初始为 1（乘法单位元）
    a %= mod;              // 先把底数取模，防止一开始就溢出
    while (n > 0) {        // 只要指数还没处理完
        if (n & 1)         // 检查 n 的最低位是不是 1
            res = res * a % mod;  // 如果是 1，把当前的 a 乘进结果
        a = a * a % mod;   // a 自乘，准备对应下一个二进制位
        n >>= 1;           // 把 n 右移一位，相当于除以 2（向下取整）
    }
    return res;
}

走一遍具体过程：算 $3^{13} \bmod 100$。

步骤	n（二进制）	n&1	res	a
初始	1101	—	1	3
第1轮	1101	1	1×3=3	3²=9
第2轮	110	0	3（不变）	9²=81
第3轮	11	1	3×81=243→43	81²→6561→61
第4轮	1	1	43×61=2623→23	61²→3721→21

最终 $3^{13} \bmod 100 = 23$。你可以验证：$3^{13}=1594323$，$1594323 \bmod 100 = 23$。✓

逆元

乘法逆元的定义：如果 $ax \equiv 1 \pmod{m}$，那 $x$ 就是 $a$ 的逆元，记作 $a^{-1}$。

逆元不一定存在。一个关键定理：$a$ 在模 $m$ 下有逆元，当且仅当 $\gcd(a, m) = 1$（即 $a$ 和 $m$ 互质）。

为什么？如果 $\gcd(a, m) = g > 1$，那 $ax$ 永远是 $g$ 的倍数，而 $ax - 1$ 不是，所以 $ax \equiv 1 \pmod{m}$ 无解。反过来，如果 $\gcd(a, m) = 1$，02-01 的扩展欧几里得告诉我们存在 $x, y$ 使得 $ax + my = 1$，对 $m$ 取模就得到 $ax \equiv 1 \pmod{m}$。

当模数 $p$ 是素数时，$\gcd(a, p) = 1$ 对所有 $a \not\equiv 0$ 都成立，所以每个非零数都有逆元——这就是为什么题目喜欢用素数模数（$10^9+7$ 或 $998244353$）。

方法一：费马小定理（模数是素数时）

费马小定理：若 $p$ 是素数且 $\gcd(a,p)=1$，则 $a^{p-1} \equiv 1 \pmod{p}$。

两边乘 $a^{-1}$：$a^{-1} \equiv a^{p-2} \pmod{p}$。

所以逆元就是一次快速幂：

long long inv(long long a, long long p) {
    return qpow(a, p - 2, p);
    // 比如算 3 在 mod 5 下的逆元：
    // qpow(3, 3, 5) = 3³ mod 5 = 27 mod 5 = 2  ✓
}

有了逆元，模意义下的除法就是乘以逆元：$\frac{b}{a} \bmod p = b \cdot a^{p-2} \bmod p$。

方法二：扩展欧几里得（模数任意时）

当模数 $m$ 不是素数时（比如 $m = 10^9+9$ 不是素数，或者 $m$ 是合数），费马小定理不能用。但只要 $\gcd(a, m) = 1$，就可以用扩展欧几里得求逆元：

// 求 ax ≡ 1 (mod m) 的 x，返回 -1 表示无解
long long mod_inv(long long a, long long m) {
    long long x, y;
    long long g = exgcd(a, m, x, y);    // ① ax + my = g
    if (g != 1) return -1;                // ② gcd ≠ 1，逆元不存在
    return (x % m + m) % m;               // ③ 保证结果在 [0, m) 内
}

• ① 扩展欧几里得求出 $x, y$ 使得 $ax + my = \gcd(a, m)$。
• ② 如果 $\gcd(a, m) \neq 1$，逆元不存在。
• ③ $ax \equiv \gcd(a,m) \pmod{m}$，当 $\gcd = 1$ 时 $ax \equiv 1$，$x$ 就是逆元。可能为负数，所以加 $m$ 再取模。

批量预处理：阶乘与逆阶乘

后面做组合计数时要频繁算 $\binom{n}{k} \bmod p$，如果每次都调快速幂就太慢了。可以预处理：

（下面以 $p = 10^9+7$ 为例；如果题目要求 $998244353$，把 MOD 改掉即可。）

const int MOD = 1e9 + 7;
const int MAXN = 3000006;
long long fact[MAXN], inv_fact[MAXN];

void init(int n) {
    fact[0] = 1;                        // 0! = 1
    for (int i = 1; i <= n; i++)
        fact[i] = fact[i-1] * i % MOD;  // fact[i] = i!

    inv_fact[n] = qpow(fact[n], MOD - 2, MOD);  // (n!)^{-1}，一次快速幂
    for (int i = n - 1; i >= 0; i--)
        inv_fact[i] = inv_fact[i+1] * (i+1) % MOD;
        // 倒推：(k!)^{-1} = ((k+1)!)^{-1} × (k+1)
        // 为什么？因为 (k+1)! = (k+1) × k!
        // 所以 k! = (k+1)! / (k+1)
        // 逆元：(k!)^{-1} = ((k+1)!)^{-1} × (k+1)
}

long long C(int n, int k) {
    if (k < 0 || k > n) return 0;
    return fact[n] * inv_fact[k] % MOD * inv_fact[n-k] % MOD;
    // C(n,k) = n! / (k! × (n-k)!)
}

重点在于倒推：不需要对每个 $k!$ 都算一次快速幂，从 $(n!)^{-1}$ 出发，乘 $(n, n-1, \ldots, 1)$ 就能得到所有逆阶乘。这样 $n$ 次快速幂变成了 1 次快速幂 + $O(n)$ 次乘法。

例题

例题 1：M&A 050 — Power（★1）

题目：给定 $a$ 和 $b$，求 $a^b \bmod (10^9+7)$。

数据范围：$1 \le a \le 100$，$1 \le b \le 10^9$

—— AtCoder M&A 050

思路：快速幂模板题。

代码：

#include <cstdio>
using namespace std;

const long long MOD = 1000000007;

int main() {
    long long a, b;
    scanf("%lld%lld", &a, &b);

    long long base = a, ans = 1;
    while (b > 0) {
        if (b & 1) ans = ans * base % MOD;
        base = base * base % MOD;
        b >>= 1;
    }
    printf("%lld\n", ans);
}

---

例题 2：Tessoku Book — A29 Power

题目：和例题 1 一样——求 $a^b \bmod (10^9+7)$。

数据范围：同 M&A 050

—— AtCoder Tessoku Book A29

（注：这道题和 M&A 050 指向同一道题，作为快速幂的第二次练习。）

---

例题 3：M&A 053 — Sum of 4^N（★2）

题目：求 $4^0 + 4^1 + \cdots + 4^N \bmod (10^9+7)$。

数据范围：$1 \le N \le 10^{18}$

—— AtCoder M&A 053

思路：等比数列求和 = $(4^{N+1} - 1) \times 3^{-1} \bmod p$。快速幂算 $4^{N+1}$，费马小定理算 $3^{-1}$。这道题在 02-03 中已经讲过。

代码：（同 02-03 例题 4，不再重复）

---

例题 4：M&A 062 — Teleporter（★3）

题目：有 $N$ 个小镇，编号 $1$ 到 $N$。小镇 $i$ 有一个传送器，传送到小镇 $A_i$。从小镇 $1$ 出发，使用传送器恰好 $K$ 次后，到达哪个小镇？

数据范围：$2 \le N \le 2 \times 10^5$，$1 \le K \le 10^{18}$

—— AtCoder M&A 062

思路：$K$ 高达 $10^{18}$，绝对不能模拟。但传送路径一定会形成环——$N$ 个小镇，从任意起点出发最多 $N$ 步就会进入循环。

核心思想：倍增（doubling）。预处理 next[i][k] 表示从小镇 $i$ 出发，使用 $2^k$ 次传送器后到达的小镇。这和快速幂的思想完全一致——把 $K$ 写成二进制，用预处理好的 $2^k$ 跳转表逐步跳。

预处理时间 $O(N \log K)$，每次查询 $O(\log K)$。

代码：

#include <cstdio>
using namespace std;

int main() {
    int N;
    long long K;
    scanf("%d%lld", &N, &K);
    int A[200010];
    for (int i = 1; i <= N; i++)
        scanf("%d", &A[i]);

    // ① 倍增表：next[i][k] = 从 i 出发 2^k 步后到达的小镇
    int next[200010][62];  // log2(10^18) ≈ 60
    for (int i = 1; i <= N; i++)
        next[i][0] = A[i];    // ② 2^0 = 1 步
    for (int k = 0; k < 60; k++)
        for (int i = 1; i <= N; i++)
            next[i][k+1] = next[next[i][k]][k];  // ③ 倍增递推

    // ④ 用快速幂的思想，把 K 分解成 2 的幂次之和
    int cur = 1;
    for (int k = 0; k < 61; k++) {
        if ((K >> k) & 1)    // ⑤ K 的第 k 位是 1
            cur = next[cur][k];  // ⑥ 跳 2^k 步
    }
    printf("%d\n", cur);
}

逐行解析：

• ② next[i][0] = A[i]：1 步（$2^0$）的转移就是直接传送。
• ③ next[i][k+1] = next[next[i][k]][k]：$2^{k+1} = 2^k + 2^k$。先跳 $2^k$ 步，再从那里跳 $2^k$ 步。这和快速幂的 base = base * base 完全一样！
• ⑤⑥ 把 $K$ 写成二进制，逐位检查。如果第 $k$ 位是 1，就跳 $2^k$ 步。最终 cur 就是 $K$ 步后的位置。
• 验证：样例 1，$N=4, K=5, A=[3,2,4,1]$。路径 $1 \to 3 \to 4 \to 1 \to 3 \to 4$，$K=5$ 后到达 4。✓

---

例题 5：M&A 055 — Recurrence Formula 1（★2）

题目：数列 $a_1 = 1, a_2 = 1, a_n = 2a_{n-1} + a_{n-2}$（$n \ge 3$）。求 $a_N \bmod (10^9+7)$。

数据范围：$3 \le N \le 10^{18}$

—— AtCoder M&A 055

思路：$N$ 最大 $10^{18}$，逐项递推 $O(N)$ 不可能。需要矩阵快速幂。

递推式 $a_n = 2a_{n-1} + a_{n-2}$ 可以写成矩阵形式：

$\begin{pmatrix} a_n \\ a_{n-1} \end{pmatrix} = \begin{pmatrix} 2 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} a_{n-1} \\ a_{n-2} \end{pmatrix}$

所以 $\begin{pmatrix} a_N \\ a_{N-1} \end{pmatrix} = M^{N-2} \begin{pmatrix} a_2 \\ a_1 \end{pmatrix}$，其中 $M$ 是上述 $2 \times 2$ 矩阵。$M^{N-2}$ 用矩阵快速幂 $O(\log N)$ 算出。

代码：

#include <cstdio>
using namespace std;

const long long MOD = 1000000007;

struct Mat {
    long long a[2][2];
    Mat() { a[0][0]=a[0][1]=a[1][0]=a[1][1]=0; }
};

Mat mul(Mat A, Mat B) {
    Mat C;
    for (int i = 0; i < 2; i++)
        for (int k = 0; k < 2; k++)
            for (int j = 0; j < 2; j++)
                C.a[i][j] = (C.a[i][j] + A.a[i][k] * B.a[k][j]) % MOD;
    return C;
}

Mat power(Mat A, long long e) {
    Mat ans;
    ans.a[0][0] = ans.a[1][1] = 1;  // ① 单位矩阵
    while (e > 0) {
        if (e & 1) ans = mul(ans, A);
        A = mul(A, A);
        e >>= 1;
    }
    return ans;
}

int main() {
    long long N;
    scanf("%lld", &N);

    Mat M;
    M.a[0][0] = 2; M.a[0][1] = 1;    // ② 递推矩阵
    M.a[1][0] = 1; M.a[1][1] = 0;

    Mat R = power(M, N - 2);          // ③ M^{N-2}

    // ④ a_N = R[0][0] * a_2 + R[0][1] * a_1
    long long ans = (R.a[0][0] * 1 + R.a[0][1] * 1) % MOD;
    printf("%lld\n", ans);
}

逐行解析：

• ② 矩阵 $M = \begin{pmatrix} 2 & 1 \\ 1 & 0 \end{pmatrix}$，满足 $M \cdot (a_{n-1}, a_{n-2})^T = (a_n, a_{n-1})^T$。
• ③ $M^{N-2}$ 用矩阵快速幂，和普通快速幂结构完全一样，只是乘法换成了矩阵乘法。
• ④ $a_N = M^{N-2}[0][0] \times a_2 + M^{N-2}[0][1] \times a_1 = R[0][0] + R[0][1]$。

验证：$N = 10$。数列为 $1, 1, 3, 7, 17, 41, 99, 239, 577, 1393$。$a_{10} = 1393$。✓

参考文献

理论讲义 — AtCoder Algorithm Lectures

• 合同式の基礎（難易度 1）
• Euclid の互除法（難易度 2）

教材讲解 — 競技プログラミングの鉄則 第 5 章

• 5.4 B29 Power Hard（快速幂 解说，含 Power 函数实现）

基础练习 — アルゴリズムと数学 演習問題集

• 050 Power（快速幂 A^B mod p）【例题】
• 053 Sum of 4^N（快速幂求等比数列和）【例题】
• 062 Teleporter（倍增/快速幂 K≤10^18）【例题】
• 055 Recurrence Formula 1（矩阵快速幂）【例题】
• 051 Combination Hard（逆元求组合数）
• 054 Fibonacci Hard（矩阵快速幂 mod 10^9）
• 056 Recurrence Formula 2（矩阵快速幂）

系统练习 — 競技プログラミングの鉄則

• A29 Power（快速幂，同M&A 050）【例题】
• A30 Combination（逆元求C(n,r) mod p）

---

系列索引

第零章 基础工具

• 00-01 复杂度与暴力枚举
• 00-02 排序
• 00-03 累积和与差分

第一章 搜索技术

• 01-01 二分搜索
• 01-02 双指针与尺取法

第二章 数学基础

• 02-01 整数与整除
• 02-02 素数与筛法
• 02-03 同余与模运算
• 02-04 快速幂与逆元
• 02-05 组合计数
• 02-06 容斥原理
• 02-07 欧拉函数
• 02-08 博弈论基础

第三章 数据结构

• 03-01 栈、队列与单调栈
• 03-02 堆与优先队列
• 03-03 并查集
• 03-04 树状数组
• 03-05 线段树
• 03-06 懒标记线段树
• 03-07 Sparse Table 与倍增
• 03-08 字符串哈希

第四章 图论

• 04-01 图的遍历
• 04-02 最短路—Dijkstra 与 01-BFS
• 04-03 最短路—Bellman-Ford 与 Floyd
• 04-04 拓扑排序
• 04-05 最小生成树
• 04-06 强连通分量与 2-SAT
• 04-07 二分图与网络流
• 04-08 树上问题

第五章 动态规划

• 05-01 DP入门—状态与转移
• 05-02 背包问题族
• 05-03 LIS、LCS与编辑距离
• 05-04 区间DP
• 05-05 状态压缩DP
• 05-06 树形DP与数位DP
• 05-07 矩阵快速幂与线性递推

第六章 贪心

• 06-01 贪心策略与交换论证
• 06-02 区间调度与优化

第七章 字符串

• 07-01 字符串哈希详述
• 07-02 后缀数组

第八章 进阶

• 08-01 卷积与FFT
• 08-02 期望与概率
• 08-03 Zeta / Möbius 变换
• 08-04 计算几何基础
• 08-05 采样与随机化