阶段 8 · 数学 · 第 42 章

快速幂与取模

★ 关键一步是 `a^b = (a^(b/2))²` —— 接第 12 章的分治,一句话把 b 次乘法压到 log b 次。⚠ 可这一章一半的内容在**取模**上:乘法溢出、负数取模、以及 `p = 1` 和 `b = 0` 这两个边界。

例题:多组 a^b mod p 建议用时:110 分钟
算法五行就写完了 —— 而这一章七成的力气花在别处

a^b mod p 的正解只有五行。可这一章要讲的东西里,算法只占三成

要讲的占多少
a^b = (a^(b/2))²(第 6 步)三成 —— 接第 12 章的分治,一句话的事
乘法溢出(p 到 10¹⁸ 时 a*b 装不下)★ 三成
负数取模(C++ 的 % 会给负余数)★ 两成
p = 1b = 0 这两个边界★★ 两成

★ 第 40 章立过一条:算法越短,风险越往「取模、溢出、负数、0 和 1」上转移。 这一章是那句话最彻底的一次现场 —— 六个错误版本里,四个和算法本身无关

⚠ 而生成器那一节会看到一个前面十五章都没出现过的形状: 有两个 bug 要「两个边界同时出现」 —— 靠各自独立的概率相乘,等一辈子也等不来。

1 一句话问题

给定 q 组询问(q ≤ 10⁵),每组三个整数 abp

  • −10¹⁸ ≤ a ≤ 10¹⁸0 ≤ b ≤ 10¹⁸1 ≤ p ≤ 10¹⁸
  • 输出 a^b mod p,答案取在 [0, p) 里;约定 0⁰ = 1

★ 最后再输出一行:有多少组的答案是 0

★ 三个范围都是故意写成这样的
范围为什么这么写卡住哪个 bug
a 可以是负数C++ 的 % 对负数给负余数wrongNeg
p 可以到 10¹⁸res * a 就是 10³⁶ —— long long 装不下wrongOverflow
p 可以等于 1这时候任何数都 ≡ 0,连 a⁰ 也是 0wrongOne
b 可以等于 0 且约定 0⁰ = 1「底数是 0 就返回 0」这句顺手的剪枝会吞掉它wrongZero

★ 而末行那句「有多少组的答案是 0」是「题面多问一句」的第八次 —— 这次它兜的是答案恰好为 0 的那些组p = 1a ≡ 0), 也就是最容易被一句「顺手的特判」搞错的那一片。

2 手算一遍:默认那 8 组

★ 这 8 组是特意排的:六个错误版本全部现形
2 10 1000000007            → 1024
-3 5 1000000007            → 999999764     ★ 负底数 + 奇次幂:(-3)⁵ = -243 ≡ p-243
5 0 1000000007             → 1             a⁰ = 1
0 0 7                      → 1             ★ 题面约定 0⁰ = 1
123456789 1000 1000000007  → 620139939
7 13 1                     → 0             ★ p = 1:任何数都 ≡ 0
999999937 12 999999999999999989 → 860406115710957185   ★ p ≈ 10¹⁸:乘法必须 __int128
5 0 1                      → 0             ★★ p = 1 且 b = 0 —— 这一组只有它抓得到 wrongOne
末行                        → 2             (答案是 0 的有两组)

⚠ 第二行值得单独看:(-3)⁵ = -243,而 C++ 里 -243 % (10⁹+7) 给的是 −243, 不是 999999764少写一句 if (a < 0) a += p;,答案就是个负数。

⚠⚠ 而最后一组是专门为 wrongOne 加的p = 1 那一行(第 6 组)它其实是对的 —— 因为 b = 13 ≥ 1,循环里 a = 7 % 1 = 0, res 乘着乘着就变 0 了。只有 p = 1b = 0 同时成立,它才现形。 ★ 这件事第 9 步会变成整整一档生成器。

3 暴力:老老实实乘 b 次

brute.cpp标准答案:连乘 b 次(和快速幂是完全不同的想法)
输入(stdin)
输出
点「运行 ▶」看结果
⚠⚠ 这一章的对拍有一个前面十五章都没碰到过的硬限制

标准答案是「连乘 b 次」,所以对拍的数据里 b 必须压得很小(这里压在 10⁴ 以内)。 可题面允许 b 到 10¹⁸ ——

★★ 那一半的取值范围,对拍根本够不着。

⇒ 补救办法在第 10 步:不比答案,比性质a^(b₁+b₂) ≡ a^b₁ · a^b₂)。 这是第 31 章「答案不唯一时写验证器」那一招的另一种用法: ★ 这一次答案是唯一的,够不着的是「规模」。

4 慢在哪:次数是能精确算出来的

★★ 这一章的次数不是「大约 log」,是三个等号
count.cpp三种做法的次数(全部是算出来的,不是跑出来的)
输出
点「运行 ▶」看结果

b = 10¹⁸(二进制 60 位,其中 24 个 1)、p ≈ 10¹⁸(二进制 60 位):

做法次数是快速幂的几倍
连乘 b 次1 000 000 000 000 000 0001.19 × 10¹⁶
快速幂 + 龟速乘(加法次数)5 04060 ×
✓ 快速幂 + __int128乘法次数)841 ×

★★ 那个 84 是这么来的,而且是个等号

乘法次数 = popcount(b)(乘进答案的次数)+ 二进制位数(平方的次数)= 24 + 60 = 84。 (和第 38 章树状数组那个「步数恰好等于 popcount(r)」是同一种数法。)

⚠ 而这张表里有一件很容易被忽略的事: 快速幂的乘法次数只和 b 有关,和 p 一点关系都没有; 只有龟速乘那一行跟着 p 走(每次乘法换成 ⌊log₂p⌋+1 次加法)。

两条曲线的旋钮根本不是同一个 —— 下一步那张耗时表会把这件事量出来。

本机实测q = 10⁵ 组询问):

g++ -O2 -o genBig genBig.cpp && g++ -O2 -o fast fast.cpp && g++ -O2 -o mul mul.cpp && g++ -O2 -o brute brute.cpp
./genBig 100000 1 > big1.txt     # 0: b≈10¹⁸ p=10⁹+7 / 1: b≈10¹⁸ p≈10¹⁸ / 2: b≤10⁴
time ./fast < big1.txt > /dev/null
数据fast__int128mul(龟速乘)brute(连乘)
b ≈ 10¹⁸,p = 10⁹+70.05 秒1.46 秒跑不动
b ≈ 10¹⁸,p ≈ 10¹⁸0.05 秒2.91 秒跑不动
b ≤ 10⁴,p = 10⁹+70.02 秒0.31 秒2.66 秒

★★ 第一行到第二行:fast 一动不动(0.05),mul 翻了一倍(1.46 → 2.91) —— 因为 p 从 30 位变成 60 位,而龟速乘的开销正比于 p 的位数。 这就是上面那句「两条曲线的旋钮不是同一个」的实测版。

genBig.cpp(三个档位)两个旋钮:b 有多大、p 有多大 —— 它们管的不是同一件事

5 ⚠ 三件和算法无关的事,先说清楚

★ ① 乘法溢出:不是「p 大」,是「p 大过一条具体的线」

resa 都在 [0, p) 里,所以 res * a 最大是 (p−1)²。 long long 的上限是 2⁶³ − 1 ≈ 9.22 × 10¹⁸,于是安全线是

p ≤ √(2⁶³) ≈ 3.04 × 10⁹

★★ 而 OI 里最顺手的模数 10⁹+7 正好在这条线以下(10⁹)² = 10¹⁸ < 9.2×10¹⁸)—— 所以拿 10⁹+7 造数据,wrongOverflow 永远抓不到。

⚠ 第 38、40 章那条:先估一估多大的数据才溢出,再决定拧到哪一格。 这一章的「那一格」精确到了小数点后两位。

两种兜法:

  • __int128(GCC 扩展,OI 里可以用)—— 正解用的就是它;
  • 龟速乘:没有 __int128 的时候,把乘法也拆成二进制。
mul.cpp⚠ 不是 bug:龟速乘版,答案逐字节相同,只是每次乘法变成 log p 次加法
输入(stdin)
输出
点「运行 ▶」看结果

★ 注意它和快速幂是同一个想法用了两遍: 快速幂把指数拆成二进制,龟速乘把乘数拆成二进制。

★ ② 负数取模:C++ 的 % 会给负余数

-7 % 3 == -1(不是 2)。C++ 规定 % 的结果和被除数同号。

⇒ 所以读进来第一件事就是

a %= p;
if (a < 0) a += p;      // 或者写成 a = ((a % p) + p) % p

两步都不能少:先 % p 是为了让它落进 (−p, p),再 + p 才能落进 [0, p)。 (只写 a += pa = −10¹⁸p = 7 时根本不够。)

★★ ③ 初值不是 1,是 1 % p

a⁰ = 1 没错,可题面要的是 a^b mod p —— 而 p = 1 时任何数都 ≡ 0。

ll res = 1 % p;      // ★ 不是 1

★ 这一句是这一章最容易漏、也最容易被小数据蒙混过去的一行: p = 1b ≥ 1 时它其实错不了(循环里 a 会变成 0,res 跟着变 0)—— 只有 p = 1b = 0 同时成立,它才现形。

6 ★ 关键一步:a^b = (a^(b/2))²

★★ 两种说法,同一件事

递归的说法(接第 12 章的分治):

  • b 是偶数:a^b = (a^(b/2))²
  • b 是奇数:a^b = a · (a^((b−1)/2))²

⇒ 每一步把指数砍一半,所以只要 ⌊log₂b⌋ + 1 步。

迭代的说法(更好记,也是代码里那五行)—— 把 b 看成二进制

b = 13 = 1101₂   ⇒   a¹³ = a⁸ · a⁴ · a¹

于是「从低位往高位扫 b,位是 1 就把当前的 a^(2^k) 乘进答案」:

ll res = 1 % p;
while (b) {
    if (b & 1) res = (i128)res * a % p;
    a = (i128)a * a % p;
    b >>= 1;
}

★ 这两种说法的次数是同一个:popcount(b) 次「乘进答案」+ 位数次「平方」

fast.cpp正解:快速幂(算法本身五行,另外三行全是取模的坑)
输入(stdin)
输出
点「运行 ▶」看结果

7 动画一:b 的二进制一位位消失

b 的二进制从低位往高位扫:位是 1 就乘进答案,然后 a 自平方
连乘要 13 次 / 快速幂 7 次
第 1 / 6 步
★ 乘进答案 0平方 0popcount(b) = 3,位数 = 4
b =
1
1
0
1
res = 1
a  = 3
★ ★ 初值是 1 % p,不是 1 —— p = 1 时答案是 0。
起点:res = 1 mod 1000000007 = 1,a = 3,b = 13(二进制 1101)。
怎么看这个动画
  • 上面那排格子是 b 的二进制(高位在左),处理过的位会灰掉 —— 这就是「每一步把指数砍一半」的画面版
  • 这一位是 1 就把当前的 a(也就是 a^(2^k))乘进 res;不管是 0 是 1,a 都要自平方。
  • ★ 右边两个计数器:乘进答案的次数 = popcount(b),平方的次数 = 位数 —— 都是等号。
  • 切到「★ b = 255(二进制全是 1)」和「★ b = 256(只有一个 1)」对比一下: 同样是 8~9 位,乘进答案的次数差了 8 倍。
  • ⚠ 再切到「p = 1」那一档:从第 0 帧起 res 就是 0 —— 初值 1 % p 那一句的画面版

trace.cpp 是这个动画的文字版,check:viz 拿它和动画逐步比。

trace.cpp动画照着它画:每一步的位、res、a,以及两个计数器
输出
点「运行 ▶」看结果

8 ★ 对拍:六个错误版本

★ 清单:两个白送,四个要「特定的值」
错误版本靠什么现形
wrongLoop while (b > 1)只要 b ≥ 1 —— ★ 基线
wrongOrder 先平方再判位同上,几乎白送
wrongNeg 不归一化负数★ 要 a < 0指数是奇数(两件事同时)
wrongOverflow 乘法没用 __int128★ 要 p > 3.04 × 10⁹
wrongOne 初值写成 1★★ 要 p = 1 且 b = 0
wrongZero0⁰ 当成 0★★ 要 a ≡ 0 且 b = 0

★★ 最后两条要的是「两个边界同时出现」—— 这是前十五章都没碰到过的形状。

wrongLoop.cpp✗ while (b > 1):最后一位没乘进去(基线)
输入(stdin)
输出
点「运行 ▶」看结果
wrongOrder.cpp✗ 先平方、再判这一位 —— 整整错一位
输入(stdin)
输出
点「运行 ▶」看结果
wrongNeg.cpp✗ 没归一化负数 —— 答案跑到 [0,p) 外面去了
输入(stdin)
输出
点「运行 ▶」看结果
wrongOverflow.cpp✗ 乘法没用 __int128 —— 只在 p > 3.04×10⁹ 时现形
输入(stdin)
输出
点「运行 ▶」看结果
wrongOne.cpp✗ 初值写成 1 —— 只在 p = 1 且 b = 0 时现形
输入(stdin)
输出
点「运行 ▶」看结果
wrongZero.cpp✗ 「底数是 0 就返回 0」—— 顺手加的一句,把 0⁰ = 1 也吞了
输入(stdin)
输出
点「运行 ▶」看结果
对拍器
★ 这个生成器的顺手档,六个 bug 里有四个是 0 / 300 —— 全书最惨的一次开局。而救活它们的办法,有一个是前面十五章都没用过的:把两个边界配着塞。

300 轮实测(种子 1..300,最终档 12):

故意写错的地方被抓第几轮
wrongLoopwhile (b > 1)278 / 300第 1 轮
wrongOrder(先平方再判位)278 / 300第 1 轮
wrongOverflow(没用 __int128278 / 300第 1 轮
wrongOne(初值写成 1)185 / 300第 1 轮
wrongNeg(不归一化负数)125 / 300第 1 轮
wrongZero0⁰ 当成 0)111 / 300第 8 轮
mul(只是慢)0 / 300

9 ★★ 生成器:四个 0,以及一个前面没用过的招

★★★ 顺手档:六个 bug 里四个 0 / 300 —— 全书最惨的一次开局
档位相对档位 0 改了什么LoopOrderNegOvfOneZero
0(顺手写法)q ∈ [3,8]、a ∈ [1,10⁹]、b ∈ [1,10⁴]、p = 10⁹+73003000000
1a 允许取负数300300244000
2★ b 允许取 03003000000
3★ a 允许取 03003000000
4p 随机取到 10¹⁸300300030000
5p 专门塞 12992990000
6p 取小([1,50])2752900000

① 那四个 0 各有各的原因,而且每一个都指着「顺手写法」的一条:

  • a 顺手取正 ⇒ wrongNeg 见不到光;
  • b 顺手从 1 起、a 顺手不取 0 ⇒ wrongZero 见不到光;
  • ★★ p 顺手写 10⁹+7(OI 里的条件反射)⇒ wrongOverflow 正好安全落地 —— 它要的不是「p 大」,是 p 大过 3.04×10⁹ 那条线,而 10⁹+7 就在线下。

② ⚠⚠ 而档位 5 是这一章最该盯的一行:专门塞了 p = 1wrongOne 还是 0 / 300。 因为它要的是 p = 1b = 0 同时成立 —— 而这一档里 b 还是从 1 起的。

单独补一个边界不够,得把两个边界配起来。

★★ 于是有了这一章的新招:把两个边界配着塞
档位内容LoopOrderNegOvfOneZero最弱支
91+2+3+4+5(五个边界全放开)287287136287664646
12(最终档)9 + ★★ 成对造边界(25% 的询问直接造成 b = 0,再配上 a ≡ 0p = 1278278125278185111111
13(对照)12 减去「成对造边界」(=档位 9)287287136287664646
10(对照)9 减去「p 塞 1」2992991812990310
11(对照)9 减去「p 到 10¹⁸」295295157062240

① 「成对造」把最弱支从 46 抬到 111(One 66 → 185、Zero 46 → 111)。

道理很简单,算一下就明白:把两个边界各自独立地塞进去, b = 0 的概率是 20%、a ≡ 0 的概率是 15% —— 同时发生只有 3%,一组询问平均 5.5 条,也就是每轮 15% 左右。 而「成对造」直接把这两件事绑在一起发生。

★★★ 有些 bug 要的不是「某一个特定的值」(第 41 章那条), 而是「两个特定的值同时出现」—— 靠各自独立的概率相乘,等一辈子也等不来。

② 两个对照档说明另外两处也不可替代: 撤掉「p 塞 1」→ wrongOne 回到 0;撤掉「p 到 10¹⁸」→ wrongOverflow 回到 0

③ ⚠ 而「成对造」是有代价的:Loop / Order / Neg / Ovf 四列各掉 10 上下 (287 → 278、136 → 125……),因为四分之一的询问被「b = 0」占掉了,那些询问测不到别的。

★ 又是第 30 章那条:把 0 变成非 0,掉多少抓获率都划算。

gen.cpp(十四个档位)七处改动全部可重跑,包括那个「单独补一个边界不够」的诊断档

10 ★★ 对拍够不着的那一半:比性质,不比答案

★ 标准答案跑不到 b = 10¹⁸,那一半只能换个验法
verify.cpp★ 验证器:b 到 10¹⁸、p 到 10¹⁸,查三条性质
输出
点「运行 ▶」看结果

查的是三条答案该满足的性质(全部用 __int128 兜住乘法):

性质式子
① 指数可加a^(b₁+b₂) ≡ a^b₁ · a^b₂ (mod p)
② 底数可乘(a₁·a₂)^b ≡ a₁^b · a₂^b (mod p)
③ 费马小定理(p 质数、a 不是 p 的倍数)a^(p−1) ≡ 1 (mod p)

两万组随机数据、b 跑到 10¹⁸ 量级(对拍那一侧只到 10⁴,差了十四个数量级)——三条全过。

⚠⚠ 但验证器的盲区必须写清楚(第 31 章那条):

这三条性质证明不了答案一定对。 比如一份「把结果统一乘 2」的程序能过 ①(不能过 ③); 而一份「永远输出 0」的程序 ①② 都能过(0 ≡ 0 · 0)——只有 ③ 抓得住它。

验证器是「够不着」时的替代品,不是对拍的升级版。

★ 而这次「够不着」的原因和第 31 章不一样,值得并排记一次:

  • 第 31 章:答案不唯一,没法逐字节比;
  • 这一章:答案是唯一的,可标准答案跑不到那个规模

11 动画二:三种做法的次数 —— 换 p 试试

三种做法的次数(★ 而它们的答案完全相同)
b 的二进制有 60 位、 其中 24 个 1; p 的二进制有 30
连乘 b 次
1,000,000,000,000,000,000
快速幂 + 龟速乘(加法次数)
2,520
✓ 快速幂 + __int128(乘法次数)
84
★★ 快速幂的次数 = popcount(b) + 位数 = 24 + 60 = 84
⚠ 条的长度按 对数 画(不然 10¹⁸ 那一根会把别的挤没);右边的数字才是真的次数。
★ 换 p 试试:**上面第一根和最后一根一动不动** —— 快速幂的乘法次数只和 b 有关;只有龟速乘那一行跟着 p 走。
★ 这个动画有一个专门要你动手的地方

b 固定、只换 p

第一根(连乘)和最后一根(快速幂)一动不动,只有中间那根(龟速乘)跟着变。

因为快速幂的乘法次数只和 b 有关 —— popcount(b) + 位数,和模数一点关系都没有; 而龟速乘把每一次乘法拆成 ⌊log₂p⌋+1 次加法,它才是跟着 p 走的那一个

⚠ 而这三根条对应的三份代码,答案逐字节相同check:viz 里 300 轮钉着)——

第 36 章那个「对拍看不见慢」的第七次现场。

12 自测

自测清单0 / 10
这一章我卡在哪(过一个月回来看,这几行比整章正文都值钱)
这一章记住三句话
  1. a^b = (a^(b/2))² —— 把 b 看成二进制,乘法次数正好是 popcount(b) + 位数。 ★ 而这个次数只和 b 有关,和模数一点关系都没有
  2. 算法五行,坑有三处:乘法溢出(p > 3.04×10⁹ 就要 __int128)、 负数取模(% 会给负余数)、初值 1 % p(p = 1 时答案是 0)。 ★ 第 40 章那句在这一章最彻底:算法越短,风险越往边界上转移。
  3. 有些 bug 要「两个边界同时出现」 —— p = 1b = 0a ≡ 0b = 0。 ⚠ 各自独立地塞进去,同时发生的概率只有 3%;必须配着造(最弱支 46 → 111)。

⚠ 下一章(组合数)会把这一章的快速幂直接拿去用:求乘法逆元。 而杨辉三角那条递推,会和第 17 章那个「记忆化搜索的调用次数正好是杨辉三角」呼应收尾。