【算法】大数乘法问题及其高效算法

题目

编写两个任意位数的大数相乘的程序，给出计算结果。好比：

题目描述： 输出两个不超过100位的大整数的乘积。
输入： 输入两个大整数，如1234567 和 123
输出： 输出乘积，如：151851741

或者

求 1234567891011121314151617181920 * 2019181716151413121110987654321 的乘积结果

分析

所谓大数相乘（Multiplication algorithm），就是指数字比较大，相乘的结果超出了基本类型的表示范围，因此这样的数不可以直接作乘法运算。

参考了不少资料，包括维基百科词条Multiplication algorithm，才知道目前大数乘法算法主要有如下几种思路：

1. 模拟小学乘法：最简单的乘法竖式手算的累加型；
2. 分治乘法：最简单的是Karatsuba乘法，通常化之后有Toom-Cook乘法；
3. 快速傅里叶变换FFT：（为了不精度问题，能够改用快速数论变换FNTT），时间复杂度O(N lgN lglgN)。具体可参照Schönhage–Strassen algorithm；
4. 中国剩余定理：把每一个数分解到一些互素的模上，而后每一个同余方程对应乘起来就行；
5. Furer’s algorithm：在渐进意义上FNTT还快的算法。不过好像不太实用，本文就不做介绍了。你们能够参考维基百科Fürer’s algorithm

解法

咱们分别实现一下以上算法，既然不能直接使用乘法作运算，最简单最容易想到的办法就是模拟乘法运算。

一、模拟乘法手算累加

 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()]);
}

看了以上的代码，感受思路虽然很简单，可是实现起来却很麻烦，那么咱们有没有别的方法来实现这个程序呢？答案是有的，接下来我来介绍第二种方法。

二、模拟乘法累加 - 改进

简单来讲，方法二就是先不算任何的进位，也就是说，将每一位相乘，相加的结果保存到同一个位置，到最后才计算进位。

例如：计算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位），最后从右向左累加时就多了一个空间。

三、分治 - Karatsuba算法

以上两种模拟乘法的手算累加型算法，他们都是模拟普通乘法的计算方式，时间复杂度都是O(n^2)，而这个Karatsuba算法，时间复杂度仅有 $O(n^{\log _{2}3})$ 。下面，我就来介绍一下这个算法。

Karatsuba于1960年发明在 $O(n^{\log _{2}3})$ 步骤内将两个n位数相乘的Karatsuba算法。它反证了安德雷·柯尔莫哥洛夫于1956年认为这个乘法须要 $\Omega (n^{2})$ 步骤的猜测。

首先来看看这个算法是怎么进行计算的，见下图：

图中显示了计算5678 * 1234的过程，首先是拆分红abcd四个部分，而后分别计算ac, bd, (a+b)*(c+d)，最后再用第三个算式的结果减去前面两个（其实获得的就是bc+ad，可是减小了乘法步骤），而后，计算式1后面加4个0，计算式2后面不加，计算式3后面加2个0，再把这三者相加，就是正确结果。

接下来，就来证实一下这个算法的正确性。这是一幅来自Karatsuba Multiplication Algorithm – Python Code的图，咱们来看看：

咱们假设要相乘的两个数是x * y。咱们能够把x，y写成：

$$x = a * 10^{n/2} + b$$

$$y = c * 10^{n/2} + d$$

这里的n是数字的位数。若是是偶数，则a和b都是n/2位的。若是n是奇数，则你可让a是n/2+1位，b是n/2位。（例如a = 12，b = 34；a = 123，b = 45），那么x*y就变成了：

$$x*y = (a * 10^{n/2}+b) * (c * 10^{n/2}+d)$$

进一步计算，

$$x*y = a * c * 10^n + (a * d + b * c) * 10^{n/2} + bd$$

对比以前的计算过程。结果已经呼之欲出了。这里惟一须要注意的两点就是：

1. (a * d + b * c)的计算为了防止两次乘法，应该使用以前的计算
2. 这些乘法在算法里应该是递归实现的，数字很大时，先拆分，而后拆分出来的数字仍是很大的话，就继续拆分，直到a * b已是一个很是简单的小问题为之。这也是分治的思想。

---

咱们举例来尝试一下这种算法，好比计算12345 * 6789，咱们让a = 12，b = 345。同时c = 6，d = 789。也就是：

$$12345 = 12 · 1000 + 345\\ 6789 = 6 · 1000 + 789$$

那么a*c，b*d的结果以下：

$$\begin{align*} & z_2 = a*c = 12 × 6 = 72\\ & z_0 = b*d = 345 × 789 = 272205 \\ & z_1 = (12 + 345) × (6 + 789) − z_2 − z_0 = 283815 − 72 − 272205 = 11538 \\ \end{align*}$$

最终结果就是：

$$\begin{align*} & result = z_2 · 10^{2*3} + z_1 · 10^3 + z_0 \\ & result = 72 · 10^6 + 11538 · 10^3 + 272205 = 83810205. \\ \end{align*}$$

以上就是使用分治的方式计算乘法的原理。上面这个算法，由 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);
}

总结：

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]

参考资料

• 维基百科：Multiplication algorithm
• 大数乘法，在算法上，主要有几种思路？
• 华为OJ机试题目：两个大整数相乘(纯C语言实现两个大整数相乘，两种方法实现大数相乘)
• 算法理解 — 大数相乘问题
• KaraTsuba乘法 — 高效的大数乘法
• Karatsuba Multiplication Algorithm – Python Code
• Toom-Cook 3-Way Multiplication from GMP Documentations
• 大数相乘(from滴答滴答百度空间)
• 大整数运算