大數乘法問題及其高效演算法
題目
編寫兩個任意位數的大數相乘的程式,給出計算結果。比如:
題目描述: 輸出兩個不超過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
解法
我們分別實現一下以上演算法,既然不能直接使用乘法做運算,最簡單最容易想到的辦法就是模擬乘法運算。
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()]); }
看了以上的程式碼,感覺思路雖然很簡單,但是實現起來卻很麻煩,那麼我們有沒有別的方法來實現這個程式呢?答案是有的,接下來我來介紹第二種方法。
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[]
陣列的最後一位空置,不可能被佔用,我們就響應地把num1的第i位與num2的第j位相乘,結果應該存放在結果的第i+j位上的這個結果往後順移一位(放到第i+j+1位
),最後從右向左累加時就多了一個空間。
3、分治 - Karatsuba演算法
Karatsuba於1960年發明將兩個n位數相乘的Karatsuba演算法。它反證了安德雷·柯爾莫哥洛夫於1956年認為這個乘法需要 $ \Omega (n^{2})$ 步驟的猜想。
首先來看看這個演算法是怎麼進行計算的,見下圖:
根據上面的思路,實現的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 就合適了,不用再去弄更復雜的遞迴乘法。
其中,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 演算法計算大數乘積。