从KaraTsuba算法谈JS的大数乘法（上）

引

之前在某社区看到JKolmogorov 和 Karatsuba 关于乘法算法的故事，遂探索一番，发现其是关于高效大数乘法的算法。众所周知，一般的乘法的时间复杂度是O(n²)，而 Karatsuba 提出的算法复杂度仅有O(N^log₂3)，社区中少有基于JS的算法实现，本文试图从Javascript语言出发谈论大数乘法。

正

分析

所谓大数乘法（Multiplication algorithm）是指数字比较大，相乘数为整数，且相乘的结果超出了基本类型的表示范围，所以这样的数不能够直接使用语言内置的乘法运算符计算得出结果。

目前大数乘法算法有以下几种思路：

1. 模拟小学乘法
2. 分治乘法：最简单的是 Karatsuba 乘法，一般化以后有 Toom-Cook 乘法
3. 快速傅里叶变换FFT：（为了避免精度问题，可以改用快速数论变换FNTT），时间复杂度O(n(logn)(loglogn))。可参照Schönhage–Strassen algorithm
4. 中国剩余定理：把每个数分解到一些互素的模上，然后每个同余方程对应乘起来就行
5. Furer’s algorithm：在渐进意义上FNTT还快的算法。

算法实现[JS]

1. 正常模拟小学乘法

乘法运算可以分拆为两步：
- 第一步，是将乘数与被乘数逐位相乘；
- 第二步，将逐位相乘得到的结果，对应相加起来。

代码如下：

function bigNumberMultiplyInMinor (a, b, baseDigit = 10) {
    const aArr = a.split('')
    const bArr = b.split('')

    let result = []

    for (let i = aArr.length - 1; i >= 0; i--) {
        let p = 0, singleArr = []

        for (let j = bArr.length - 1; j >= 0; j--) {
            const r = aArr[i] * bArr[j] + p

            singleArr.unshift(r % baseDigit)
            p = Math.floor(r / baseDigit)
        }

        if (p !== 0) singleArr.unshift(p)

        if (i === aArr.length - 1) {
            result = singleArr.slice()
        } else {
            let carry = 0
            let idx = result.length - 1 - (aArr.length - 1 - i)

            while(!!singleArr.length) {
                const singleVal = singleArr.pop()

                if (idx < 0) {
                    result.unshift(singleVal + carry)
                    carry = 0
                } else {
                    const resultVal = result[idx]
                    const val = singleVal + resultVal + carry
                    const newVal = val % baseDigit

                    carry = Math.floor(val / baseDigit)
                    result[idx] = newVal
                }

                idx--
            }

            if (carry !== 0) {
                result.unshift(carry)
            }
        }
    }

    return result.join('')
}

2. 模拟乘法【改进】

我们分析乘法过程，发现两个数相乘的结果的长度总不能超过相乘数累加的长度。
举例：“1234(length = 4) x 56(length = 2) = 69104(length = 5 < 4 + 2)”

            1234
            * 56
   -------------
   0  6 12 18 24
   5 10 15 20  0
   -------------
   5 16 27 38 24
   -------------
   6  9  1  0  4

所以我们先抛弃固有的小学乘法的思维模式，将两数按位相乘，并对同一位置的结果进行累加，最后处理进位，也可得到正确的值。
也就是说，将每一位相乘，相加的结果保存到同一个位置，到最后才计算进位。
代码如下：

function bigNumberMultiplyInMinorImprove (a, b) {
    const aArr = a.split('')
    const bArr = b.split('')
    const aLen = aArr.length
    const bLen = bArr.length

    let result = new Array(aLen + bLen).fill(0)

    for (let i = 0; i < aLen; i++) {
        for (let j = 0; j < bLen; j++) {
            result[i+j+1] += (aArr[i] * bArr[j]) 
        }
    }

    for (let k = result.length - 1; k >0; k--) {
        if (result[k] >= 10) {
            result[k-1] += Math.floor(result[k] / 10)
            result[k] %= 10
        }
    }

    return result.join('').replace(/^0*(?=\d)/g, '')
}

可以看到这个版本的算法实现很简洁。当然，算法的时间复杂度还是O(n²)。

3. 分治 - Karatsuba 算法

所以我们说了这么长的时间，Karatsuba 算法到底是什么呢？
上面两种模拟乘法都是手算累加型算法，时间复杂度都是O(n²)，Karatsuba算法是一种快速乘法算法。 它是由 Anatoly Karatsuba 于1960年发现并于1962年发表的。它将两个n位数的乘法减少到最多n^log₂⁡3≈n1.585一般的单位数乘法。 因此，它比传统算法更快，后者需要n²个单位数乘积。Karatsuba 算法使用了分治的思想来解决这个问题。

我们假设要相乘的两个数是x * y。我们可以把x，y写成：

x = a ∗ 10 ^n/2 + b y = c ∗ 10 ^n/2 + d

我们定义halfN的值为上面的n/2，n是数字的位数，如果是偶数，则分割后到的低位数和高位数都是n/2位的，如果n是奇数，则让高位数是n/2+1位，低位数是n/2位。那么x*y则可表示为：

x * y = (a * 10 ^n/2 + b) * (c * 10 ^n/2 + d) 进而得到： x * y = a * c * 10 ⁿ + (a * d + b * c) * 10 ^n/2 + bd

而 (a * d + b * c) 的结果可由以下公式得出：

(a + b) * (c + d) - a * c - b * d 上面就是分治的乘法原理，上面这个算法，由 Anatolii Alexeevitch Karatsuba 于1960年提出并于1962年发表，所以也被称为 Karatsuba 乘法。

这些乘法在算法里应该是递归实现的，数字很大时，先拆分，然后拆分出来的数字还是很大的话，就继续拆分，直到a * b已经是一个非常简单的小问题为之。这也是分治的思想。

代码实现如下：

function karatsuba(strA, strB) {
    if (Number(strA) === 0 || Number(strB) === 0) return 0
    if (Number(strA) <= 10 && Number(strB) <= 10) return strA * strB

    const aLen = strA.length
    const bLen = strB.length
    const halfN = Math.floor(Math.max(aLen, bLen) / 2)

    const firstA = aLen > halfN ? strA.slice(0, aLen - halfN) : '0'
    const lastA = aLen > halfN ? strA.slice(aLen - halfN) : strA
    const firstB = bLen > halfN ? strB.slice(0, bLen - halfN) : '0'
    const lastB = bLen > halfN ? strB.slice(bLen - halfN) : strB

    // 很多位的运算会导致数值溢出，需要自己实现加法运算
    const plusA = String((Number(firstA) + Number(lastA))) 
    const plusB = String((Number(firstB) + Number(lastB)))

    const c2 = karatsuba(firstA, firstB)
    const c0 = karatsuba(lastA, lastB)
    const c1 = karatsuba(plusA, plusB) - c0 - c2
    const result = c2 * Math.pow(10, (2 * halfN)) + c1 * Math.pow(10, halfN) + c0

    return result
}

在实现这个算法过程中，遇到几个点需要说一下：

1. halfN的取值问题，在Java8的源码中，halfN是这样定义的：int half = (Math.max(xlen, ylen)+1)/2，我的疑惑在于如果两个乘数相差位数很大，比如说：12345678910 * 110，这个时候，halfN的值为5,超过了第二位乘数的位数，如果处理？只需让c的值为0，然后d的值是原始的数。以上面的举例来说，c为0，d为110.
2. 在分割乘数时，我们需保证分割后的乘数都是同样的位移数。其次，对于分割后的低位数和高位数的位数不同的情况，需要保证位数最长的那个乘数分割后的高位数的位数比低位数要多。举例说明：有这样两个乘数：12345 * 110，根据上面的定义，n为5，halfN则为2，所以第一个乘数的分割结果是：高位数是123，低位数是45，第二个乘数的分割结果是：高位数是1，低位数是10。这样才能保证最后得到正确的结果。注意看分割高位数和低位数时，对slice方法的参数的处理。len - halfN就是为了让高位数保留更多的位数。
3. 可以看到，在只实现大数乘法运算中，会用到原生的运算符。这样导致在对两个很多位的数相乘时，原生的加法，乘法运算会由于数值溢出从而导致运算结果错误，所以要在Javascript中使用 Karatsuba 乘法，需实现相应的大数加法操作。

链

• 【算法】大数乘法问题及其高效算法
• Kolmogorov和Karatsuba关于乘法算法的故事
• 力扣-43.字符串相乘
• Schönhage–Strassen algorithm
• Karatsuba算法
• GNU MP 6.1.2: Karatsuba Multiplication
• 大数乘法，在算法上，主要有几种思路？ - 知乎