算法与数据结构 / 分治 · 从递归式到随机化 / 计数 · 取幂 · 大整数乘法 待审核 2 / 3
Karatsuba · 快速幂

计数 · 取幂 · 大整数乘法

分治最出名的实例是排序,但排序恰好是最不需要动脑的那一类:切成两半、各自递归、线性合并,f(n)f(n) 摆在明面上。本页的三个例子各自动了别的地方:一个在合并时多算了一样东西,一个折半的对象根本不是问题规模,一个把 aa 从 4 压到了 3。

1 · 归并顺带数出的 inversion

序列 AA 里满足 i<ji < jAi>AjA_i > A_j 的下标对叫一个 inversion。它计量的是「离排好序还有多远」:升序序列有 0 个,倒序序列有 n(n1)/2n(n-1)/2 个。二重循环数一遍是 Θ(n2)\Theta(n^2)

按下标落在哪一半,inversion 分成三类:两个下标都在左半、都在右半、一左一右。前两类就是两个规模减半的同类子问题,第三类必须在合并处算掉。

合并时两半已各自有序。从右半取走 RjR_j 的那一刻,左半还没被取走的元素全都比 RjR_j 大,而它们的原始下标又全在 RjR_j 之前,于是这一步一次记下 Li|L| - i 个 inversion。左半剩几个是现成的读数,不必再比较。整个算法的形状与 merge sort 完全一样,T(n)=2T(n/2)+Θ(n)T(n) = 2T(n/2) + \Theta(n),代价 Θ(nlogn)\Theta(n \log n)

注 · 递归回来时两半已被排过序,元素的位置全乱了,跨半区的 inversion 数却不受影响:左半内部怎么重排,都不改变左半任一元素与右半任一元素的大小关系,而跨半区 inversion 只由这种关系决定。左右两半各自内部的 inversion 已经在递归里数过,排序把它们清零并不丢账。

图 1-1 · 归并求 inversion 计数的单步执行。取右半元素时左半的剩余长度一次性计入,可换样例观察升序与倒序两个极端。

警示 · 比较写成 LiRjL_i \le R_j 走左半,不能写成 <<。相等的两个元素不构成 inversion,<< 会在遇到相等时改走右半,把左半剩余长度错记一笔;同一处改动还会让归并失去稳定性。classics.test.ts[2, 2, 2, 2] 应得 0 这条断言守的就是它。

2 · 折半的是指数

ana^n 连乘要 n1n-1 次乘法。折半的对象换成指数:an=(an/2)2a^n = (a^{n/2})^2,奇数时多乘一个 aa。递归式是

T(n)=T(n/2)+Θ(1)T(n) = T(n/2) + \Theta(1)

对照 分治的骨架与复杂度 §3.2 的表,a=1a = 1b=2b = 2c=0c = 0:分水岭 nlog21=n0=1n^{\log_2 1} = n^0 = 1ff 同阶,属情形二,答案 Θ(logn)\Theta(\log n)。每层只剩一个子问题,log\log 因子全部来自层数。

自底向上的写法把递归摊平:按指数的二进制位从低到高走,维护一个不断自乘的 base,遇到为 1 的位就把它乘进累加器。乘法次数等于二进制位里 1 的个数加上平方次数,也就是 log2n\lfloor \log_2 n \rfloor 加上 popcount。powFramesn=2301n = 2^{30} - 1 上用掉 59 次乘法,朴素连乘要 1073741822 次。

这套折半只依赖乘法的结合律,与被乘的东西是什么无关。换成 2×22 \times 2 矩阵,就得到线性递推的对数解法:

(1110)n=(F(n+1)F(n)F(n)F(n1))\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n = \begin{pmatrix} F(n+1) & F(n) \\ F(n) & F(n-1) \end{pmatrix}

F(40)F(40) 只需 7 次矩阵乘法,而逐项递推要走 40 步。差距在 nn 大时才真正显出来:n=1018n = 10^{18} 的 Fibonacci 取模,线性递推走不完,矩阵快速幂是 60 次矩阵乘法。

图 2-1 · 快速幂按指数的二进制位单步执行。可在标量与 Fibonacci 转移矩阵之间切换,两侧的乘法次数完全一致。

警示 · 标量版的模数取的是 1000003,不是竞赛里常见的 1000000007。JS 的 number 是 double,精确整数只到 2539.007×10152^{53} \approx 9.007 \times 10^{15};两个接近 10910^9 的数相乘会到 101810^{18},超出后不报错,只是静默丢掉低位。取 106+310^6 + 3 之后乘积不到 101210^{12},仍在精确区间内。要用大模数就得换 BigInt 或把乘法拆成高低位两段。

3 · Karatsuba 的第三次乘法

两个 nn 位十进制数相乘,竖式要做 n2n^2 次一位数乘法。把两个数各从中间切开,记 B=10n/2B = 10^{n/2}

x=xhB+xl,y=yhB+ylx = x_h B + x_l, \quad y = y_h B + y_l
xy=xhyhB2+(xhyl+xlyh)B+xlylxy = x_h y_h B^2 + (x_h y_l + x_l y_h) B + x_l y_l

直接照这个式子算要四次半长乘法,T(n)=4T(n/2)+Θ(n)T(n) = 4T(n/2) + \Theta(n)。分水岭是 log24=2\log_2 4 = 2,比 c=1c = 1 大,属情形一,答案仍是 Θ(n2)\Theta(n^2)。切了等于白切。

Karatsuba 的观察是中间那一项不必单独乘出来。展开

(xhxl)(ylyh)=xhyl+xlyhxhyhxlyl(x_h - x_l)(y_l - y_h) = x_h y_l + x_l y_h - x_h y_h - x_l y_l

右边后两项正是已经算出的 xhyhx_h y_hxlylx_l y_l,把它们加回去就得到中间项。三次半长乘法加上若干次加减法:

T(n)=3T(n/2)+Θ(n)T(n) = 3T(n/2) + \Theta(n)

分水岭降到 log231.585\log_2 3 \approx 1.585,仍属情形一,答案变成 Θ(n1.585)\Theta(n^{1.585})ff 一动没动,省下来的全部来自 aa 从 4 变成 3。

karatsuba 实测的一位数乘法次数:4 位 9 次对竖式的 16 次,8 位 27 次对 64 次,16 位 81 次对 256 次。倍率恰好是每加倍一次位宽,乘法次数乘 3 而非乘 4。

图 3-1 · Karatsuba 递归树的先序单步执行。每个内部结点分出三个半宽子问题,可换被乘数观察叶子数始终是 3 的位宽对数次幂。

注 · 教科书多写成 (xh+xl)(yh+yl)(x_h + x_l)(y_h + y_l) 的加法变体,实现里换成了上面的减法变体,原因是递归树的形状。两个 n/2n/2 位数相加会进位到 n/2+1n/2 + 1 位,那一支递归下去每层都要多带一位,树就不再齐整,叶子数也对不上 3log2n3^{\log_2 n}。差值 xhxlx_h - x_l 则始终落在 (B,B)(-B, B) 内,位宽严格减半,代价是要处理负数的符号。classics.test.ts 直接断言了这个结构:8 位输入下叶子恰好 27 个、内部结点 13 个、树高 3。

警示 · 预设里的 8 位被乘数是挑过的。整个实现走的是 double,而 8 位数最大可到 99999999,平方后约 101610^{16},已经越过 2532^{53}。预设取 12345678 与 87654321,乘积约 1.08×10151.08 \times 10^{15},仍然精确。要跑真正的大整数得把数字换成数组表示,那时基例的「一位数乘法」才当得起这个名字。

Strassen 的矩阵乘法是同一个套路换了个舞台:把 n×nn \times n 矩阵切成四个 n/2n/2 阶块,朴素做法要八次块乘(log28=3\log_2 8 = 3,答案 n3n^3),Strassen 用七次加上大量加减法(log272.807\log_2 7 \approx 2.807)。压的同样是 aa

4 · 参考文献

  1. Karatsuba, A., & Ofman, Y. (1962). Multiplication of many-digital numbers by automatic computers. Doklady Akademii Nauk SSSR, 145(2), 293–294.
  2. Strassen, V. (1969). Gaussian elimination is not optimal. Numerische Mathematik, 13(4), 354–356.
  3. Knuth, D. E. (1997). The Art of Computer Programming, Volume 2: Seminumerical Algorithms (3rd ed.). Addison-Wesley.
  4. Kleinberg, J., & Tardos, É. (2005). Algorithm Design. Addison-Wesley.