大数乘法(一)|大数乘法中的分治法

636 阅读3分钟

之前在快速幂中讲过,复杂的快速幂时间复杂度会受到大数乘法的限制,这次我是来填坑的!我才不会告诉你之前我本来只想讲个分治结果被卷到了才找的奇技淫巧快速幂,嗯,这次就满足我的心愿:讲一次分治嘿嘿~

010.png

大数乘法一览

大数乘法包括

  • 小学模拟乘法:最简单的乘法竖式手算累加型;O(n2)O(n^2)
  • Schönhage–Strassen 史恩哈格·施特拉森算法,复杂度为 Θ(nlog(n)log(logn))\Theta(n \log( n) \log (\log n))
  • Karatsuba (卡拉楚巴)快速乘法算法 (大数乘法)(分治最简单的)Onlog23O(n^{\log_23})
  • Toom-Cook乘法/Toom-3算法 (比Karatsuba更有效的大数乘法)(也是分治)O(nlog35)O(n^{\log_3 5})
  • 快速傅里叶变换FFT(不是分治)Onlogn)O(n\log n) 渐进最优 asymptotically optimal

Integer multiplication in time O(nlogn)O(n \log n) David Harvey and Joris van der Hoeven
高斯重采样技术,将整数乘法问题转化为一系列在复数域上的多维离散傅里叶变换,其维数都是2的幂。这些变换可以通过Nussbaumer的快速多项式变换快速求解。

  • 富尔算法

  • 中国剩余定理:把每个数分解到一些互素的模上,然后每个同余方程对应乘起来就行;

etc.

本次讲解大数乘法中的分治法: Karatsuba && Toom-Cook (卡拉楚巴和图姆库克)

为啥要写上汉字音译!当然是为了防止装逼失败啊!(并没有)

Karatsuba ( Toom-Cook 在 2 的特例 )

概述

Karatsuba算法主要应用于两个大数的相乘,原理是将大数分成两段后变成较小的数位,然后做3次乘法,并附带少量的加法操作和移位操作。 

步骤

现有两个大数,xyx,y。 

首先将xyx,y分别拆开成为两部分,可得x1x0y1y0x_1,x_0,y_1,y_0。他们的关系如下: 

x=x1×10m+x0x = x_1 \times 10m + x_0; 

y=y1×10m+y0y = y_1 \times 10m + y_0。其中m为正整数,m<nm < n,且x0y0x_0,y_0 小于 10m。 

那么 

xy=(x1×10m+x0)(y1×10m+y0)xy = (x_1 \times 10^m + x_0)(y_1 \times 10^m + y_0)

=z2×102m+z1×10m+z0\quad=z_2 \times 10^{2m} + z_1 \times 10^m + z_0,其中: 

        z2=x1×y1z_2 = x_1 \times y_1

        z1=x1×y0+x0×y1z_1 = x_1 \times y_0 + x_0 \times y_1; 

        z0=x0×y0z_0 = x_0 \times y_0。 

此步骤共需4次乘法,但是由Karatsuba改进以后仅需要3次乘法。因为: 

z1=x1×y0+x0×y1 z_1 = x_1 \times y_0+ x_0 \times y_1 

z1=(x1+x0)×(y1+y0)x1×y1x0×y0z_1 = (x_1 + x_0) \times (y_1 + y_0) - x_1 \times y_1 - x_0 \times y_0,

x0×y0x_0 \times y_0 便可以由加减法得到。 

所以: 

x×y=z2×102m+z1×10m+z0x\times y = z_2 \times 10^{2m} + z_1 \times 10^m + z_0 

z2=x1×y1 z_2 = x_1 \times y_1 

z1=(x1+x0)×(y1+y0)x1×y1x0×y0=(x1+x0)×(y1+y0)x1×y1z0z_1 = (x_1 + x_0) \times (y_1 + y_0) - x_1 \times y_1 - x_0 \times y_0 \\ \quad \quad \quad = (x1 + x_0) \times (y_1 + y_0) - x_1 \times y_1 - z_0

z0=x0×y0 z_0 = x_0 \times y_0 

Recursively computer 递归计算(x1×y1)(x_1\times y_1) 

Recursively computer 递归计算(x1+x0)×(y1+y0) (x_1 + x_0) \times (y_1 + y_0) 

Recursively computer 递归计算(x0×y0) (x_0 \times y_0) 

实例讲解

x=12345y=6789x = 12345,y=6789,令m=3m=3。那么有: 

12345=12×1000+34512345 = 12 \times 1000 + 345 

6789=6×1000+7896789 = 6 \times 1000 + 789 

下面计算: 

z2=12×6=72z_2 = 12 \times 6 = 72 

z0=345×789=272205z_0 = 345 \times 789 = 272205

z1=(12+345)×(6+789)z2z0=11538z_1 = (12 + 345) \times (6 + 789) - z_2 - z_0 = 11538

然后我们按照移位公式xy=z2×102m+z1×101m+z0(xy = z_2 \times 10^2m + z_1 \times 10^1m + z_0)可得: 

xy=72×10002+11538×1000+272205=83810205xy = 72 \times 10002 + 11538 \times 1000 + 272205 = 83810205

Toom-Cook (又称 Toom-3)

概述

图姆-库克算法的原理是:

对于给定的两个大整数 aa bb, 将  aa bb 分成 kk个较小的部分, 每个部分的长度为 ll, 并对这些部分执行运算。随着 kk 的增长, 可以组合许多乘法子运算, 从而降低算法的整体复杂度, 然后再次使用图姆-库克算法递归计算乘法子运算, 依此类推。 

Toom-3和图姆-库克两个术语有时会被错误的 混用, 但事实上Toom-3只是图姆-库克算法在 k=3k=3时的特例。

Toom-3将9次乘法降低至仅需5次, 使其在 Θ(nlog(5)/log(3))Θ(n1.46)\Theta\left(n^{\log (5) / \log (3)}\right) \approx \Theta\left(n^{1.46}\right) 的时间里运行。 

步骤

假设我们希望将 A=123,456,789A=123,456,789B=987,654,321B=987,654,321相乘。每个整数被分成三个等长的数字,ll位宽。     

l=3l=3 时,我们可以用两个多项式表示两个操作数: 

 a(x)=123x2+456x+789  a(x) = 123x^2 + 456x + 789 

 b(x)=987x2+654x+321  b(x) = 987x^2 + 654x + 321 

其中 A=a(1,000)A=a(1,000)B=b(1,000)B=b(1,000)。 

为了计算 c(x)=a(x)×b(x)c(x) = a(x) \times b(x)

我们在五个点,例如 {2,1,0,1,2}\{-2, -1, 0, 1, 2\}上评估两个多项式 aabb,并将结果相乘

问题:为什么是5个点?
因为x2x^2 相乘,我们知c(x)c(x) 必为4阶,所以选(n+1=)5(n+1 = )5个点
有限差分法,仅使用加法即可计算出 A 和 B 的连续值(但这里不是) 

x=2时,{A=3692369B=296122961AB=1092609  \quad x= -2时, \\ \left\{\begin{array}{l} A= 369 & (-2 , 369)\\ B= 2961 &(-2,2961) \\ AB= 1092609 \end{array} \right.

x=1时,{A=4561456B=6541654AB=298224   \quad x= -1时, \\ \left\{\begin{array}{l} A= 456 & (-1 , 456)\\ B= 654 &(-1,654) \\ AB= 298224  \end{array} \right.

x=0时,{A=7890789B=3210321AB=253269   \quad x= 0时, \\ \left\{\begin{array}{l} A= 789 & (0 , 789)\\ B= 321 &(0,321) \\ AB= 253269  \end{array} \right.

x=1时,{A=136811368B=196211962AB=2684016  \quad x= 1时, \\ \left\{\begin{array}{l} A= 1368 & (1 , 1368)\\ B= 1962 &(1,1962) \\ AB= 2684016 \end{array} \right.

x=2时,{A=219322193B=557725577AB=12230361   \quad x= 2时, \\ \left\{\begin{array}{l} A= 2193 & (2 , 2193)\\ B= 5577 &(2,5577) \\ AB= 12230361  \end{array} \right.     

根据五个点的值找到多项式 P(x)=A(x)B(x)P(x)=A(x)B(x) 

P(x)=p4x4+p3x3+p2x2+p1x+p0P(x) = p_4x^4 + p_3x^3 + p_2x^2 + p_1x + p_0

我们可以代入这些点处的 x 值,得到一组在未知数 pip_i 中呈线性的联立方程: 

{16p48p3+4p22p1+p0=109260916p4p3+p2p1+p0=298224p0=253269p4+p3+p2+p1+p0=268401616p4+8p3+4p2+2p1+p0=1223036116 \left\{\begin{array}{r} 16p_4 - 8p_3 + 4p_2 - 2p_1 + p_0 = &109260916\\ p_4 - p_3 + p_2 - p_1 + p_0 = &298224\\p_0 = &253269\\p_4 + p_3 + p_2 + p_1 + p_0 = &2684016\\16p_4 + 8p_3 + 4p_2 + 2p_1 + p_0 = &1223036116 \end{array}\right.

列为线性方程,解得  {p4=121401p3=530514p2=1116450p1=662382p0=253269\left\{\begin{array}{l}p_4 = 121401\\p_3 = 530514\\ p_2 = 1116450\\ p_1 = 662382\\p_0 = 253269 \end{array}\right.

x=1000x=1000代入P(x)=p4x4+p3x3+p2x2+p1x+p0P(x) = p_4x^4 + p_3x^3 + p_2x^2 + p_1x + p_0中,

也即每个移位三位后相加得到

P(1000)=AB=121,932,631,112,635,269P(1000) = A \cdot B = 121,932,631,112,635,269

 

时间复杂度分析

通常, Toom- kk 的时间复杂度为 Θ(c(k)ne)\Theta\left(c(k) n^e\right),

其中

e=log(2k1)/log(k)e=\log (2 k-1) / \log (k) 。 

nen^e 是在乘法子运算上花费的时间,

cc 则是花费在对小常数进行的加法和乘法运算上的时间。 

大数乘法的对比综述

著名的Karatsuba算法实际上是图姆-库克算法的特例, 在Karatsuba算法中, 原始乘数被拆分成两个较小的数, 而原本的4次乘法运算缩减为 3 次, 使之在 Θ(nlog(3)/log(2))Θ(n1.58)\Theta\left(n^{\log (3) / \log (2)}\right) \approx \Theta\left(n^{1.58}\right)的时间 内完成运算。 

Toom-1等价于普通的长乘法, 具有Θ(n2)\Theta\left(n^2\right) 的复杂度。 

尽管可以通过增加 k 来使指数 e 任意接近1, 但函数增长速度非常快。 

混合级别图姆-库克算法的增长率直到 2005 年仍然是一个广为研究的开放性问题 

根据高德纳所描述算法的一种实现, 其复杂度可降低至 Θ(n22lognlogn)\Theta\left(n^{ 2^{\sqrt{2 \log n}}} \log n\right)

由于工作时的开销, 当乘数包括较小的数时, 图姆-库克算法会比长乘法更慢, 因此它适用于中等规模的乘法。 

经验结果表明,精心设计的 Toom 算法将在除最庞大的计算之外的所有方面击败 FFT 方法。规模增加会导致收益递减,虽然该方法可以接近具有较大 kk 值的快速傅立叶变换方法的 O(nlogn)O(n log n),但乘法之外的计算工作会迅速变得非常重要。 

对于更大规模的数据, 则有渐进更快的史恩哈格·施特拉森算法,复杂度为 Θ(nlognloglogn)\Theta(n \log n \log \log n)

参考资料