Karastuba算法实现高精度乘法

博主之前的文章 [C#]高精度算法 中，简述过基于循环高精度加法来实现高精度乘法的思路。
这一篇，将从另一个角度来实现这个功能。

背景 ＆ 思想

Karastuba 博士在 1960 年提出了一个用于计算大整数乘法的算法，其思想是把两个大整数的乘法转化为若干次小规模的乘法和少量的加法。

一般来说，对于两个 n 位的大整数 x 和 y ，可以将其分解为如下的两部分：

$$x={x_1}^1{0}^{\frac{n}{2}}+x_0,0\le x_0<10^{\frac{n}{2}}$$$$y={y_1}^1{0}^{\frac{n}{2}}+y_0,0\le y_0<10^{\frac{n}{2}}$$

例如： 12345 = 123 x $10^2$ + 45
其实质是利用欧几里得算式将一个整数分解成 m = qn + r 的形式。
因此，
$x\ast y=(x_1^1{0}^{\frac{n}{2}}+x_0)(y_1^1{0}^{\frac{n}{2}}+y_0)$
$={x_1{y}_1}^1{0}^n+(x_1{y}_0+x_0{y}_1{)}^1{0}^{\frac{n}{2}}+x_0{y}_0$
这其中， $10^n$ 可以通过位移来高效运算，原先的大整数乘法就可以分解为四次小规模的乘法和一些加法。
但是，实际运算中给出的两个数字不一定是等位数的，例如： 123 x 12345 。
假设 x 和 y 分别是 m 位和 n 位的大整数，
则：
$$x=x_1^1{0}^{\frac{m}{2}}+x_0,0\le x_0<10^{\frac{m}{2}}$$$$y=y_1^1{0}^{\frac{n}{2}}+y_0,0\le y_0<10^{\frac{n}{2}}$$
易得：
$$xy=(x_{1}10^{\frac{m}{2}}+x_0)(y_{1}10^{\frac{n}{2}}+y_0)$$$$=x_1{y}_{1}10^{\frac{m+n}{2}}+x_1{y}_{0}10^{\frac{m}{2}}+x_0{y}_{1}10^{\frac{n}{2}}+x_0{y}_0$$
那么我们只要反复迭代多项式的每一项 $x_i{y}_j$ , 直到其中一个乘数只有一个为止。

实现

分析：x=123， y=321。用上标表示迭代的次数，则 x · y 的过程如下：
$123=12x10^{1}3321=32x10^{1}1$
$x_1^{①}=12,x_0^{①}=3y_1^{①}=32,y_0^{①}=1$
那么，
$xy=x_1^{①}{y}_1^{①}10^2+x_1^{①}{y}_0^{①}10^1+x_0^{①}{y}_1^{①}10^1+x_0^{①}{y}_0^{①}$
$=(12\ast 32)10^2+(12\ast 1)10^1+(3\ast 32)10^1+3\ast 1$
其中，12x1、3x32、3x1 都出现了有一个乘数是 1 的情况，那么就不必继续迭代，可以直接返回结果。
第二次迭代，实际就是对 12 、32 再次进行 karastuba 分解而已，原理是一样的。
结合思路与迭代经验，可以写出这样的代码：

def karastuba(x,y):
    if x==0 or y==0:return 0
    m,n=len(str(x)),len(str(y))
    if m==1 or n==1:return x * y
    m//=2;n//=2
    x1,x0,y1,y0=x//(10**m),x%(10**m),y//(10**n),y%(10**n)
    x1y1,x1y0,x0y1,x0y0=karastuba(x1,y1),karastuba(x1,y0),karastuba(x0,y1),karastuba(x0,y0)
    return x1y1*(10**(m+n)) + x1y0*(10**m) + x0y1*(10**n) + x0y0

使用时，调用 karastuba 函数即可。
受限于篇幅关系，支持符号和小数点的算法可以自行探索。
另外，我们可以用这样的代码来测试：

def test(x,y):
    rt = karastuba(x,y)
    print(str(x)+" * "+str(y)+" = "+str(rt), rt == x * y)
    
if __name__ == "karastuba":
    test(123,321)
    test(123456789,987654321)

最后，保存，运行，输出结果：

>>> from karastuba import *
123 * 321 = 39483 True
123456789 * 987654321 = 121932631112635269 True

参考资料

《程序员数学从零开始》，孙博 (@我是 8 位的)