跳至主要內容

Karatsuba大整数乘法

周子力大约 6 分钟教学文档算法设计与分析分治法递归大整数乘法高难度Python

Karatsuba大整数乘法

题目描述

给定两个任意精度的整数 x、y,使用 Karatsuba 分治算法求其精确乘积。允许正数、负数和0,非整数输入抛出 TypeError。

要求将较大整数拆成高位、低位,递归时用三次较小整数乘法重建乘积。为保证边界清晰,本实现约定两个非负操作数均不超过8个二进制位时,直接使用 Python 的整数乘法;更大的乘法必须按 Karatsuba 公式分解。

示例: 1234 × 5678 返回7006652,-1234 × 5678 返回-7006652,任意整数乘0返回0。

快速幂通过复用相同半次幂减少子调用;本题通过代数恒等式,把原本需要的四个乘法子问题变为三个。这里的规模是整数的二进制位数,不是整数数值大小。

题目分析

从四次乘法到三次乘法

对非负整数,取 B=2hB=2^h,写成:

x=aB+b,y=cB+d,0≤b,d<B. x=aB+b,\qquad y=cB+d,\qquad 0\leq b,d<B.

直接展开得到:

xy=acB2+(ad+bc)B+bd. xy=acB^2+(ad+bc)B+bd.

若分别求 ac,ad,bc,bdac,ad,bc,bd,需要四次子乘法。Karatsuba 改为只求:

z2=ac,z0=bd,z1=(a+b)(c+d). z_2=ac,\qquad z_0=bd,\qquad z_1=(a+b)(c+d).

因为 z1−z2−z0=ad+bcz_1-z_2-z_0=ad+bc,所以:

xy=z2B2+(z1−z2−z0)B+z0. xy=z_2B^2+(z_1-z_2-z_0)B+z_0.

多出的加减和移位均能按位数线性完成,而较大乘法的递归数量减少了一个。

分割、边界与符号

内部函数只处理非负数。令 ww 为两数二进制位数的最大值,取 h=⌊w/2⌋h=\lfloor w/2\rfloor。右移 hh 位取得高位,与掩码 B−1B-1 按位与取得低位;乘以 BB、B2B^2 则用左移实现。

当某个操作数为0时立即返回0;当 w≤8w\leq8 时直接乘法作为常数规模边界。更大的调用中,高低部分的位数至多约一半;相加可能多一位,但 w>8w>8 时子问题仍严格变小,递归必定终止。

外层记录结果符号,对绝对值进行递归,最后恢复符号。两数同号得到非负结果,异号得到非正结果。

正确性依据

边界乘法直接正确。假设三次子乘法均正确,代入上述恒等式即可得到当前两数的精确乘积;位移只是乘以2的幂,并未舍入。最后恢复符号,也符合整数乘法的符号规则。

可结合 MIT 的 Karatsuba 课程资料open in new window学习这种减少子乘法的思路。

实现代码

Python 二进制分割实现

def karatsuba(x, y):
    if type(x) is not int or type(y) is not int:
        raise TypeError("操作数必须是整数")
    negative = (x < 0) != (y < 0)

    def multiply(a, b):
        if a == 0 or b == 0:
            return 0
        width = max(a.bit_length(), b.bit_length())
        if width <= 8:
            return a * b

        half = width // 2
        mask = (1 << half) - 1
        a_high, a_low = a >> half, a & mask
        b_high, b_low = b >> half, b & mask

        z2 = multiply(a_high, b_high)
        z0 = multiply(a_low, b_low)
        z1 = multiply(a_high + a_low, b_high + b_low)
        middle = z1 - z2 - z0
        return (z2 << (2 * half)) + (middle << half) + z0

    result = multiply(abs(x), abs(y))
    return -result if negative else result


for x, y in [(1234, 5678), (-1234, 5678), (-12, -34), (0, 5678), (1, 999)]:
    print(f"{x} × {y} = {karatsuba(x, y)}")

输出:

1234 × 5678 = 7006652
-1234 × 5678 = -7006652
-12 × -34 = 408
0 × 5678 = 0
1 × 999 = 999

代码没有用十进制字符串确定位数,而是用 bit_length();低位提取和重组也全部采用二进制操作。阈值8是为了让课堂例子能展示递归,实际工程实现通常会根据运行环境选择更大的阈值。

算法执行过程

主示例中1234有11个二进制位,5678有13个位。因此 width=13、half=6、B=64B=64:

1234=19×64+18,5678=88×64+46. 1234=19\times64+18,\qquad5678=88\times64+46.

子计算调用输入返回结果
高位乘法 z2z_2(19, 88)1672
低位乘法 z0z_0(18, 46)828
和的乘法 z1z_1(37, 134)4958
中间项4958−1672−8284958-1672-8282458
重建乘积1672×4096+2458×64+8281672\times4096+2458\times64+8287006652

三个子调用的最大操作数都不超过8位,直接触发边界乘法。根调用只做三次子乘法,没有单独计算 19×4619\times46 和 18×8818\times88。

输入为 (-1234,5678) 时,内部过程完全相同,最终给7006652加上负号。规模更大的整数则会在三条分支继续采用同样的分割规则。

复杂度分析

设两输入绝对值的最大二进制位数为 ww,令 d=max⁡(1,w)d=\max(1,w)。采用按位运算成本模型:长度为 dd 的整数加减、移位、掩码操作按 O(d)O(d) 计,固定8位以内的乘法按常数计。

时间复杂度: 最坏为 O(dlog⁡23)≈O(d1.585)O(d^{\log_2 3})\approx O(d^{1.585})。每次执行三个约半长的子乘法,当前层额外工作为线性;高低位之和可能产生进位,因此更准确的规模上界为:

T(d)≤3T(⌈d/2⌉+1)+O(d). T(d)\leq3T(\lceil d/2\rceil+1)+O(d).

这个固定的一位增量不改变渐近指数。常见的 3T(d/2)+O(d)3T(d/2)+O(d) 写法是忽略进位后的简化模型。两数位数很不均衡或含0时,实际工作可以更少。传统按位竖式乘法的最坏时间为 O(d2)O(d^2)。

辅助空间复杂度: O(d)O(d) 位的同时存活的中间整数,递归栈包含 O(log⁡(d+1))O(\log(d+1)) 个调用帧;沿一条调用链的操作数长度按几何级数递减。结果本身需要 O(d)O(d) 位。若只数整数对象而不数其位数,会低估大整数存储。

Python 内置大整数乘法自身已有优化,加上本实现的递归开销,不能仅凭渐近式断言这段教学代码会比内置 x * y 更快。

易错点

  • 仍然递归计算四个交叉乘积。 会回到四分支结构,失去 Karatsuba 减少子乘法的关键。
  • 把 z1z_1 当作中间系数。 它包含 acac 和 bdbd,必须再减去 z2,z0z_2,z_0。
  • 高位只左移一次。 z2z_2 对应 B2B^2,应左移 2 * half 位。
  • 直接对负数拆位。 负数的移位和掩码含符号语义,本实现先取绝对值再处理。
  • 用数值大小作为递推规模。 分割的是表示长度,应使用位数 dd 分析。
  • 把大整数加法视为常数,同时又讨论按位乘法。 分析必须使用一致的计数模型。

课堂练习

  1. 取十进制 B=10B=10,将12与34分割,用三次乘法得到结果。
  2. 主示例的中间系数2458,是否等于 19×46+18×8819\times46+18\times88?
  3. 同样减半,四次子乘法对应的简化递推是什么?与三次相比结论如何?
  4. 输入 (-7,0)、(-7,-9) 及 (7.0,9),分别有什么结果?
参考答案
  1. a=1,b=2,c=3,d=4a=1,b=2,c=3,d=4。z2=3,z0=8,z1=3×7=21z_2=3,z_0=8,z_1=3\times7=21,中间系数为10,结果为 3×100+10×10+8=4083\times100+10\times10+8=408。这是手算公式;代码采用二进制且此例直接触发8位边界。
  2. 是。19×46=87419\times46=874,18×88=158418\times88=1584,合计2458。
  3. T(d)=4T(d/2)+O(d)T(d)=4T(d/2)+O(d),得到 O(d2)O(d^2);改为三个子乘法得到 O(dlog⁡23)O(d^{\log_2 3})。
  4. 分别返回0、63,第三个抛出 TypeError,因为7.0为浮点数。
上次编辑于:
贡献者: zlzhou