【算法】大数乘法问题及其高效算法
一、题目
编写两个任意位数的大数相乘的程序，给出计算结果。比如：
题目描述：输出两个不超过
100位的大整数的乘积。
输入：输入两个大整数，如1234567和123
输出：输出乘积，如：151851741
或者
求1234567891011121314151617181920 * 2019181716151413121110987654321的乘积结果
二、分析
所谓大数相乘(Multiplication algorithm)，就是指数字比较大，相乘的结果超出了基本类型的表示范围，所以这样的数不能够直接做乘法运算。
参考了很多资料，包括维基百科词条Multiplication algorithm，才知道目前大数乘法算法主要有以下几种思路：
- 模拟小学乘法：最简单的乘法竖式手算的累加型；
- 分治乘法：最简单的是
Karatsuba乘法，一般化以后有Toom-Cook乘法； - 快速傅里叶变换FFT：（为了避免精度问题，可以改用快速数论变换FNTT），时间复杂度\(O(N lgN lglgN)\)。具体可参照Schönhage–Strassen algorithm；
- 中国剩余定理：把每个数分解到一些互素的模上，然后每个同余方程对应乘起来就行；
- Furer’s algorithm：在渐进意义上
FNTT还快的算法，不过好像不太实用。大家可以参考维基百科Fürer’s algorithm
三、解法
我们分别实现一下以上算法，既然不能直接使用乘法做运算，最简单最容易想到的办法就是模拟乘法运算。
3.1 模拟乘法手算累加
7 8 9 6 5 2
× 3 2 1 1
-----------------
7 8 9 6 5 2 <---- 第1趟
7 8 9 6 5 2 <---- 第2趟
.......... <---- 第n趟
-----------------
? ? ? ? ? ? ? ? <---- 最后的值用另一个数组表示
如上所示，乘法运算可以分拆为两步：
- 第一步，是将乘数与被乘数逐位相乘；
- 第二步，将逐位相乘得到的结果，对应相加起来。
这有点类似小学数学中，计算乘法时通常采用的“竖式运算”。用Java简单实现了这个算法，代码如下：
/**
* 大数相乘 - 模拟乘法手算累加
*/
public static Integer[] bigNumberMultiply(int[] arr1， int[] arr2) {
ArrayList<Integer> result = new ArrayList<>(); //中间求和的结果
//arr2 逐位与arr1相乘
for (int i = arr2.length - 1; i >= 0; i--) {
int carry = 0;
ArrayList<Integer> singleList = new ArrayList<>();
//arr2 逐位单次乘法的结果
for (int j = arr1.length - 1; j >= 0; j--) {
int r = arr2[i] * arr1[j] + carry;
int digit = r % 10;
carry = r / 10;
singleList.add(digit);
}
if (carry != 0) {
singleList.add(carry);
}
int resultCarry = 0, count = 0;
int k = 0;
int l = 0;
int offset = arr2.length - 1 - i; //加法的偏移位
ArrayList<Integer> middleResult = new ArrayList<>();
//arr2每位乘法的结果与上一轮的求和结果相加，从右向左做加法并进位
while (k < singleList.size() || l < result.size()) {
int kv = 0, lv = 0;
if (k < singleList.size() && count >= offset) {
kv = singleList.get(k++);
}
if (l < result.size()) {
lv = result.get(l++);
}
int sum = resultCarry + kv + lv;
middleResult.add(sum % 10); //相加结果从右向左（高位到低位）暂时存储，最后需要逆向输出
resultCarry = sum / 10;
count++;
}
if (resultCarry != 0) {
middleResult.add(resultCarry);
}
result.clear();
result = middleResult;
}
Collections.reverse(result); //逆向输出结果
return result.toArray(new Integer[result.size()]);
}
看了以上的代码，感觉思路虽然很简单，但是实现起来却很麻烦，那么我们有没有别的方法来实现这个程序呢？答案是有的，接下来我来介绍第二种方法。
3.2 模拟乘法累加 - 改进
简单来说，方法二就是先不算任何的进位，也就是说，将每一位相乘，相加的结果保存到同一个位置，到最后才计算进位。
例如：计算98×21,步骤如下
9 8
× 2 1
-------------
(9)(8) <---- 第1趟: 98×1的每一位结果
(18)(16) <---- 第2趟: 98×2的每一位结果
-------------
(18)(25)(8) <---- 这里就是相对位的和，还没有累加进位
这里唯一要注意的便是进位问题，我们可以先不考虑进位，当所有位对应相加，产生结果之后，再考虑。从右向左依次累加，如果该位的数字大于10，那么我们用取余运算，在该位上只保留取余后的个位数，而将十位数进位（通过模运算得到）累加到高位便可，循环直到累加完毕。
核心代码如下：
/**
* 大数相乘方法二
*/
public static int[] bigNumberMultiply2(int[] num1, int[] num2) {
// 分配一个空间，用来存储运算的结果，num1长的数 * num2长的数，结果不会超过num1+num2长
int[] result = new int[num1.length + num2.length];
// 先不考虑进位问题，根据竖式的乘法运算，num1的第i位与num2的第j位相乘，结果应该存放在结果的第i+j位上
for (int i = 0; i < num1.length; i++) {
for (int j = 0; j < num2.length; j++) {
result[i + j + 1] += num1[i] * num2[j]; // 因为进位的问题，最终放置到第i+j+1位
}
}
//单独处理进位
for (int k = result.length-1; k > 0; k--) {
if (result[k] > 10) {
result[k - 1] += result[k] / 10;
result[k] %= 10;
}
}
return result;
}
！！注意：这里的进位有个大坑，因为
result[]数组是从左到右记录相对位的和（还没有进位），而最后的进位是从右向左累加进位，这样的话，如果最高位，也就是最左侧那一位的累加结果需要进位的话，result[]数组就没有空间存放了。
而正好result[]数组的最后一位空置，不可能被占用，我们就响应地把num1的第i位与num2的第j位相乘，结果应该存放在结果的第i+j位上的这个结果往后顺移一位（放到第i+j+1位），最后从右向左累加时就多了一个空间。
3.3 分治 - Karatsuba算法
以上两种模拟乘法的手算累加型算法，他们都是模拟普通乘法的计算方式，时间复杂度都是\(O(n^2)\)，而这个Karatsuba算法，时间复杂度仅有\(O(n^{log}2^3)\)。下面，我就来介绍一下这个算法。
Karatsuba于1960年发明在\(O(n^{log}2^3)\)步骤内将两个n位数相乘的Karatsuba算法。它反证了安德雷·柯尔莫哥洛夫于1956年认为这个乘法需要\(Ω(n^2)\)步骤的猜想。
首先来看看这个算法是怎么进行计算的，见下图：
Karatsuba Multiplication Algorithm步骤
图中显示了计算5678 * 1234的过程，首先是拆分成abcd四个部分，然后分别计算ac, bd, (a+b)*(c+d)，最后再用第三个算式的结果减去前面两个（其实得到的就是bc+ad，但是减少了乘法步骤），然后，计算式1后面加4个0，计算式2后面不加，计算式3后面加2个0，再把这三者相加，就是正确结果。
接下来，就来证明一下这个算法的正确性。这是一幅来自Karatsuba Multiplication Algorithm – Python Code的图，我们来看看：
Karatsuba算法证明
我们假设要相乘的两个数是x * y。我们可以把x，y写成：
这里的n是数字的位数。如果是偶数，则a和b都是n/2位的。如果n是奇数，则你可以让a是n/2+1位，b是n/2位。（例如a = 12，b = 34；a = 123，b = 45），那么x*y就可以换算为：
注意最后一步，这个式子倒数第二步中的(ad + bc)，没必要另外进行两次乘法，可以使用((a+b)*(c+d) - ac - bd)来重复利用前面的两次乘积结果ac和bd。
也就是说，我们发现最终只需要计算三次乘法a*c，b*d，(a+b)*(c+d)以及六次加法。因此这样复杂度就变为
对比之前的计算过程，结果已经呼之欲出了。这里唯一需要注意的两点就是：
(a*d + b*c)的计算为了防止两次乘法，应该使用之前的计算也就是第一幅图第四步的③-②-①- 这些乘法在算法里应该是递归实现的，数字很大时，先拆分，然后拆分出来的数字还是很大的话，就继续拆分，直到
a * b已经是一个非常简单的小问题为之。这也是分治的思想。
我们举例来尝试一下这种算法，比如计算12345 * 6789，我们让a = 12，b = 345。同时c = 6，d = 789。也就是：
那么a*c，b*d的结果如下：
最终结果就是：
以上就是使用分治的方式计算乘法的原理。上面这个算法，由Anatolii Alexeevitch Karatsuba于1960年提出并于1962年发表，所以也被称为Karatsuba乘法。
根据上面的思路，实现的Karatsuba乘法代码如下：
/**
* Karatsuba乘法
*/
public static long karatsuba(long num1, long num2) {
//递归终止条件
if (num1 < 10 || num2 < 10) return num1 * num2;
// 计算拆分长度
int size1 = String.valueOf(num1).length();
int size2 = String.valueOf(num2).length();
int halfN = Math.max(size1, size2) / 2;
/* 拆分为a, b, c, d */
long a = Long.valueOf(String.valueOf(num1).substring(0, size1 - halfN));
long b = Long.valueOf(String.valueOf(num1).substring(size1 - halfN));
long c = Long.valueOf(String.valueOf(num2).substring(0, size2 - halfN));
long d = Long.valueOf(String.valueOf(num2).substring(size2 - halfN));
// 计算z2, z0, z1, 此处的乘法使用递归
long z2 = karatsuba(a, c);
long z0 = karatsuba(b, d);
long z1 = karatsuba((a + b), (c + d)) - z0 - z2;
return (long)(z2 * Math.pow(10, (2 * halfN)) + z1 * Math.pow(10, halfN) + z0);
}
3.4 小结
Karatsuba算法是比较简单的递归乘法，把输入拆分成2部分，不过对于更大的数，可以把输入拆分成3部分甚至4部分。拆分为3部分时，可以使用下面的Toom-Cook 3-way 乘法，复杂度降低到\(O(n^{1.465})\)。拆分为4部分时，使用Toom-Cook 4-way 乘法，复杂度进一步下降到\(O(n^{1.404})\)。对于更大的数字，可以拆成100段，使用快速傅里叶变换FFT，复杂度接近线性，大约是\(O(n^{1.149})\)。可以看出，分割越大，时间复杂度就越低，但是所要计算的中间项以及合并最终结果的过程就会越复杂，开销会增加，因此分割点上升，对于公钥加密，暂时用不到太大的整数，所以使用Karatsuba就合适了，不用再去弄更复杂的递归乘法。
四、测试程序
public class LeetCodeTest {
public static void main(String[] args) {
//String a = "1234567891011121314151617181920";
//String b = "2019181716151413121110987654321";
//String a = "999999999999";
//String b = "999999999999";
//String a = "24566";
//String b = "452053";
String a = "98";
String b = "21";
char[] charArr1 = a.trim().toCharArray();
char[] charArr2 = b.trim().toCharArray();
// 字符数组转换为int[]数组
int[] arr1 = new int[charArr1.length];
int[] arr2 = new int[charArr2.length];
for (int i = 0; i < charArr1.length; i++) {
arr1[i] = charArr1[i] - '0';
}
for (int i = 0; i < charArr2.length; i++) {
arr2[i] = charArr2[i] - '0';
}
// 开始计算
int[] result = LeetcodeTest.bigNumberMultiply2(arr1, arr2);
System.out.println(a + " * " + b + " = " + Arrays.toString(result).replace(", ", ""));
}
}
最后，是测试用例输出结果：
1234567891011121314151617181920 * 2019181716151413121110987654321
= [02492816912877266687794240983772975935013386905490061131076320]
999999999999 * 999999999999 = [999999999998000000000001]
24566 * 452053 = [11105133998]
98 * 21 = [2058]
五、Java中BigInteger的乘法实现
目前最著名的高精度整数运算库是GNU的GMP，GMP是The GNU MP Bignum Library，是一个开源的数学运算库，它可以用于任意精度的数学运算，包括有符号整数、有理数和浮点数。它本身并没有精度限制，只取决于机器的硬件情况。许多著名的计算机代数系统如Axiom,Maple,Mathematica,Maxima等的底层高精度整数运算都是基于GMP实现的。
设n为乘数的位数, 就目前已知的情况而言, 不同乘法算法的时间复杂度可以从平凡的\(O(n^2)\)(普通乘法), \(O(n^{log_{2^3}})\)(Karatsuba乘法)， \(O(n^{log_{3^5}})\) (Toom-3乘法)，\(O(nlog^∗n)\) (复数域上的FFT)，其中\(log^∗n=logn(loglogn)(logloglogn)···,\)和\(O(n(logn)(loglogn))\)(有限域上的FFT),其中“有限域上的FFT”的时间复杂度已经相当接近线性了。
但是这些乘法算法中复杂度较低的算法往往有较大的常数因子，因此如果乘数的位数较少,普通乘法反而是最快的，所以实用中常常将这些不同的乘法算法结合起来使用，每次做乘法时都根据相乘两数的大小动态地选择具体采用哪一种算法，而每种算法的最佳适用范围往往依赖于具体实现和硬件环境，因此一般直接通过实验来确定。
5.1 Java中的实现
- 在
Java 7里面，就是用二重循环直接乘的。源代码：[BigInteger - Java7]（见1165行） - 在
Java 8里面，根据两个因数的大小，有三种乘法。源代码：[BigInteger - Java8]（见1464行）：- 当两个因数均小于\(2^{32×80}\)时，用二重循环直接相乘，复杂度为\(O(n^2)\)，n为因数位数（下同）；
- 否则，当两个因数均小于\(2^{32×240}\)时，采用
Karatsuba algorithm，其复杂度为\(O(n^{log_{2^3}})≈O(n^{1.585})\)； - 否则，采用
Toom-Cook multiplication，其复杂度为\(O(n^{log^{_35}})≈O(n^{1.465})\)。
其中，Java8中的源代码如下：
private static final int MULTIPLY_SQUARE_THRESHOLD = 20;
private static final int KARATSUBA_THRESHOLD = 80;
private static final int TOOM_COOK_THRESHOLD = 240;
public BigInteger multiply(BigInteger val) {
if (val.signum == 0 || signum == 0)
return ZERO;
int xlen = mag.length;
if (val == this && xlen > MULTIPLY_SQUARE_THRESHOLD) {
return square();
}
int ylen = val.mag.length;
if ((xlen < KARATSUBA_THRESHOLD) || (ylen < KARATSUBA_THRESHOLD)) {
int resultSign = signum == val.signum ? 1 : -1;
if (val.mag.length == 1) {
return multiplyByInt(mag, val.mag[0], resultSign);
}
if (mag.length == 1) {
return multiplyByInt(val.mag, mag[0], resultSign);
}
int[] result = multiplyToLen(mag, xlen, val.mag, ylen, null);
result = trustedStripLeadingZeroInts(result);
return new BigInteger(result, resultSign);
} else {
if ((xlen < TOOM_COOK_THRESHOLD) && (ylen < TOOM_COOK_THRESHOLD)) {
// 采用 Karatsuba algorithm 算法
return multiplyKaratsuba(this, val);
} else {
// 采用 Toom-Cook multiplication 3路乘法
return multiplyToomCook3(this, val);
}
}
}
我们可以看到，Java8依据两个因数的量级分别使用Karatsuba algorithm和Toom-Cook multiplication算法计算大数乘积。
5.2 Karatsuba algorithm
/**
* Java8中的 Karatsuba algorithm 算法
*/
private static BigInteger multiplyKaratsuba(BigInteger x, BigInteger y) {
int xlen = x.mag.length;
int ylen = y.mag.length;
// The number of ints in each half of the number.
int half = (Math.max(xlen, ylen) + 1) / 2;
// xl and yl are the lower halves of x and y respectively,
// xh and yh are the upper halves.
BigInteger xl = x.getLower(half);
BigInteger xh = x.getUpper(half);
BigInteger yl = y.getLower(half);
BigInteger yh = y.getUpper(half);
BigInteger p1 = xh.multiply(yh); // p1 = xh * yh
BigInteger p2 = xl.multiply(yl); // p2 = xl * yl
// p3 = (xh+xl) * (yh+yl)
BigInteger p3 = xh.add(xl).multiply(yh.add(yl));
// result = p1 * 2^(32 * 2 * half) + (p3 - p1 - p2) * 2^(32 * half) + p2
BigInteger result = p1.shiftLeft(32 * half).add(p3.subtract(p1).subtract(p2))
.shiftLeft(32 * half).add(p2);
if (x.signum != y.signum) {
return result.negate();
} else {
return result;
}
}
5.3 Toom-Cook multiplication
/**
* Java8中的 Toom-Cook multiplication 3路乘法
*/
private static BigInteger multiplyToomCook3(BigInteger a, BigInteger b) {
int alen = a.mag.length;
int blen = b.mag.length;
int largest = Math.max(alen, blen);
// k is the size (in ints) of the lower-order slices.
int k = (largest + 2) / 3; // Equal to ceil(largest/3)
// r is the size (in ints) of the highest-order slice.
int r = largest - 2 * k;
// Obtain slices of the numbers. a2 and b2 are the most significant
// bits of the numbers a and b, and a0 and b0 the least significant.
BigInteger a0, a1, a2, b0, b1, b2;
a2 = a.getToomSlice(k, r, 0, largest);
a1 = a.getToomSlice(k, r, 1, largest);
a0 = a.getToomSlice(k, r, 2, largest);
b2 = b.getToomSlice(k, r, 0, largest);
b1 = b.getToomSlice(k, r, 1, largest);
b0 = b.getToomSlice(k, r, 2, largest);
BigInteger v0, v1, v2, vm1, vinf, t1, t2, tm1, da1, db1;
v0 = a0.multiply(b0);
da1 = a2.add(a0);
db1 = b2.add(b0);
vm1 = da1.subtract(a1).multiply(db1.subtract(b1));
da1 = da1.add(a1);
db1 = db1.add(b1);
v1 = da1.multiply(db1);
v2 = da1.add(a2).shiftLeft(1).subtract(a0)
.multiply(db1.add(b2).shiftLeft(1).subtract(b0));
vinf = a2.multiply(b2);
// The algorithm requires two divisions by 2 and one by 3.
// All divisions are known to be exact, that is, they do not produce
// remainders, and all results are positive. The divisions by 2 are
// implemented as right shifts which are relatively efficient, leaving
// only an exact division by 3, which is done by a specialized
// linear-time algorithm.
t2 = v2.subtract(vm1).exactDivideBy3();
tm1 = v1.subtract(vm1).shiftRight(1);
t1 = v1.subtract(v0);
t2 = t2.subtract(t1).shiftRight(1);
t1 = t1.subtract(tm1).subtract(vinf);
t2 = t2.subtract(vinf.shiftLeft(1));
tm1 = tm1.subtract(t2);
// Number of bits to shift left.
int ss = k * 32;
BigInteger result = vinf.shiftLeft(ss).add(t2).shiftLeft(ss)
.add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);
if (a.signum != b.signum) {
return result.negate();
} else {
return result;
}
}
