大数运算 -- 乘法

异想时刻

在做数字运算的时候，常常会遇到结果溢出的情形，也就是超出了整型所能表达的最大数字。此时便不能通过整型去表达数字了，需要以字符串的形式来存储和计算数字。以下记录下几种大数乘法的运算方法以便后续应用。

普通算法

采用按位相乘错位相加的方式计算结果。若针对两个n位数的整数，其时间复杂度为 O(n_2) 。以56*34为例。

1
2
3
4
5
6
7
8
9
     5  6
  x  3  4
---------
    20 24   //按位相乘
 15 18
---------       
 15 38 24   //错位累加
---------
1 9  0  4   //进位

Karatsuba算法

卡拉楚巴算法是一种快速乘法算法，它由阿纳托利·阿列克谢耶维奇·卡拉楚巴于1960年提出并于1962年发表。若针对两个n位数的整数，其时间复杂度为 O(3n^{log_{2}3}) 。相对于普通的经典算法，在性能上有了很大的提高。

原理剖析

以 x*y 为例，x的位数 n_1 ，y的位数 n_2 ，x和y可以表达为

x=x_0*10^m+x_1,

y=y_0*10^m+y_1,

m=n/2, n=max\{n_1, n_2\}

=>

z=x*y=x_0*y_0*10^{2m}+(x_0*y_1+y_0*x_1)*10^m+x_1*y_1

=>

z_0=x_0*y_0

z_1=x_0*y_1+x_1*y_0=(x_0+x_1)(y_0+y_1)-x_0*y_0-x_1*y_1

z_2=x_1*y_1

可以看出，只需要做三次乘法便可求出系数 z_0,z_1,z_2 。递归以上过程，直到 n_1=1||n_2=1

伪代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
procedure karatsuba(x, y) {
  if(x<10||y<10) return x*y;
  //拆分数字
  int n1 = len(x), n2=len(y);
  int n = max(n1, n2);
  int m = floor(n/2);
  x0, x1 = split(x, m);
  y0, y1 = split(y, m);
  //计算系数
  z0 = karatsuba(x0, y0);
  z2 = karatsuba(x1, y1);
  z1 = karatsuba(x0+x1, y0+y1)-z0-z2;
  return z0*10^(2*m)+z1*10^m+z2;
}

Schönhage–Strassen 算法

颂哈吉-施特拉森算法是渐近快速的大整数乘法算法，由阿诺德·颂哈吉和沃尔克·施特拉森在1971年发明。算法采用快速傅里叶变换，对数字进行数论转换，可以快速计算出两个大数的乘积。若针对两个n位数的整数，其时间复杂度为 O(n\cdot logn\cdot loglogn) 。

傅里叶变换

傅里叶变换是一种线性积分变换，常用于信号在时域和频域间的变换。傅里叶变换是可逆的，可通过转换后的复函数还原成实函数。

连续傅里叶变换

将一组函数 f 映射到另一组函数 F ，其表现形式如下

F(\xi) = a\int_{-\infty}^{\infty}f(x)e^{-2\pi ix\xi}dx, \text{x为任意实数}, a \in \{1, \frac{1}{\sqrt{2\pi}}\}

a为常规化因子，由使用者自主选择，不同领域取值不一。

逆变换形式为：

f(x)=a\int_{-\infty}^{\infty}F(\xi)e^{2\pi ix\xi}d\xi,\text{x为任意实数}, a\in \{1, \frac{1}{\sqrt{2\pi}}\}

离散傅里叶变换DFT

对于N点序列 {x[n]}_{0<=n<N} 的离散傅里叶变换为

X_k = \sum^{N-1}_{n=0}e^{-i2\pi nk/N}x_n,\text{k=0,1,...,N-1}

逆变换，

x_n=\frac{1}{N}\sum_{k=0}^{N-1}e^{2\pi ink/N}X_k,\text{n=0,1,...,N-1}

原理剖析

此算法采用了数论DFT来计算两个数的乘积。注意并不是直接使用DFT，在计算时需要进行一次数论变换。为什么？由于快速傅里叶变换是复数运算，这样就存在大量的浮点数运算。不仅会导致计算性能低，而且运算结果也会存在较大的偏差。

X·Y=IDFT(DFT(X)·DFT(Y))

经过数论转换后的DFT形式如下。

X_k = \sum^{N-1}_{n=0}e^{-i2\pi nk/N}x_n\ (mod\ M),\text{k=0,1,...,N-1}

逆变换，

x_n=\frac{1}{N}\sum_{k=0}^{N-1}e^{2\pi ink/N}X_k\ (mod\ M),\text{n=0,1,...,N-1}

其中，

1，M是一个质数

2，N是M-1的因子

3， N^{-1} \equiv N\ (mod\ M)

4， \exists \alpha, \alpha^N=1\ (mod\ M), and\ \alpha^k\neq 1\ (mod\ M), k=1,...,N-1

计算步骤

Step 1, 将相乘的两项转化为两个向量，比如1234*5678，变换后为 X=[4,3,2,1,0,0,0,0] 和 Y=[8,7,6,5,0,0,0,0] ，向量的维数为两个数的位数之和。

Step 2, 对X和Y进行快速傅里叶变换。以下为具体的转换流程。

X_0=x_0+x_1+...+x_{N-1}\ (mod\ M)

X_1=x_0+e^{-2\pi i\frac{1}{N}}x_1+e^{-2\pi i\frac{2}{N}}x_2+...+e^{-2\pi i\frac{N-1}{N}}x_{N-1}\ (mod\ M)

...

X_{N-1}=x_0+e^{-2\pi i\frac{N-1}{N}}x_1+e^{-2\pi i\frac{2(N-1)}{N}}x_2+...+e^{-2\pi i\frac{(N-1)(N-1)}{N}}x_{N-1}\ (mod\ M)

=>

\begin{bmatrix} X_0\\X_1\\...\\X_{N-1} \end{bmatrix} = \begin{bmatrix} 1&1&1&...&1\\1&w&w^2&...&w^{N-1}\\...&...&...&...&...\\1&w^{N-1}&w^{2(N-1)}&...&w^{(N-1)(N-1)} \end{bmatrix} \begin{bmatrix} x_0\\x_1\\...\\x_{N-1} \end{bmatrix}\ (mod\ M),w=e^{-2\pi i/N}

Step 3, 对X和Y按为相乘，并对M取模，得到向量Z。然后求取Z的DFT逆变换，得到向量z。

\begin{bmatrix} Z_0\\Z_1\\...\\Z_{N-1} \end{bmatrix}=\begin{bmatrix} X_0*Y_0\\X_1*Y_1\\...\\X_{N-1}*Y_{N-1} \end{bmatrix}\ (mod\ M)

z_0=\frac{1}{N}(Z_0+Z_1+...+Z_{N-1})

z_1=\frac{1}{N}(Z_0+e^{2\pi i\frac{1}{N}}Z_1+...+e^{2\pi i\frac{N-1}{N}}Z_{N-1})

...

z_{N-1}=\frac{1}{N}(Z_0+e^{2\pi i\frac{N-1}{N}}Z_1+...+e^{2\pi i\frac{(N-1)(N-1)}{N}}Z_{N-1})

=>

\begin{bmatrix} z_0\\z_1\\...\\z_{N-1} \end{bmatrix}=\frac{1}{N}\begin{bmatrix} 1&1&1&...&1\\1&w^{-1}&w^{-2}&...&w^{-(N-1)}\\...&...&...&...&...\\1&w^{-(N-1)}&w^{-2(N-1)}&...&w^{-(N-1)(N-1)} \end{bmatrix} \begin{bmatrix} Z_0\\Z_1\\...\\Z_{N-1} \end{bmatrix}\ (mod\ M),w=e^{-2\pi i/N}

Step 4, 逐项遍历z，对每一项超过10的部分进行进位累加到下一项中，直到遍历结束。将z倒置便得到最终结果。

Pic1展示了算法的整个计算过程。

Pic1, 算法计算过程

References:
[1]https://zh.wikipedia.org/zh-cn/%E5%8D%A1%E6%8B%89%E6%A5%9A%E5%B7%B4%E7%AE%97%E6%B3%95
[2]https://zh.wikipedia.org/zh-cn/%E9%A0%8C%E5%93%88%E5%90%89-%E6%96%BD%E7%89%B9%E6%8B%89%E6%A3%AE%E6%BC%94%E7%AE%97%E6%B3%95
[3]https://zh.wikipedia.org/zh-cn/%E6%95%B8%E8%AB%96%E8%BD%89%E6%8F%9B
[4]https://zh.wikipedia.org/zh-cn/%E6%A8%A1%E5%8F%8D%E5%85%83%E7%B4%A0
[5]https://zh.wikipedia.org/zh-cn/%E5%BF%AB%E9%80%9F%E5%82%85%E9%87%8C%E5%8F%B6%E5%8F%98%E6%8D%A2
[6]https://zh.wikipedia.org/zh-cn/%E7%A6%BB%E6%95%A3%E5%82%85%E9%87%8C%E5%8F%B6%E5%8F%98%E6%8D%A2
[7]https://zh.wikipedia.org/wiki/%E5%82%85%E9%87%8C%E5%8F%B6%E5%8F%98%E6%8D%A2
[8]https://zh.wikipedia.org/wiki/%E8%BF%9E%E7%BB%AD%E5%82%85%E9%87%8C%E5%8F%B6%E5%8F%98%E6%8D%A2

算法
#算法，大数乘法，big number multiplication #karatsuba #Schönhage–Strassen
大数运算 -- 乘法
https://r-future.github.io/post/big-number-multiply/
Author
Future
Posted on
July 24, 2022
Licensed under