问题描述:给定两个巨大的整数,求他们的乘积。
这个问题有个比较简单的方法来计算,即利用快速傅里叶变换。下面简述下原理。
假设我们要计算\(x*y\)的值,首先把数字表示成为多项式形式:
\(x=a_0*10^0+a_1*10^1+a_2*10^2+...+a_n*10^n\)
\(y=b_0*10^0+b_1*10^1+b_2*10^2+...+b_m*10^m\)
比如数字1234可以表示成\(4*10^0+3*10^1+2*10^2+1*10^3\)
这样,计算\(x*y\)其实就是计算两个多项式的乘积,即计算多项式各项(\(10^n\))系数。
我们知道傅里叶变换的公式如下:
\(h(x)=\int f(\tau)g(x-\tau)d\tau\)
其离散形式就是把积分变成求和。
把上述\(x*y\)按照多项式乘积计算展开可以惊喜的发现,\(x*y\)的多项式系数刚好是\(x\)的系数与\(y\)的系数的卷积。其实也很好理解,因为卷积运算本身就是移动+求和,和乘法公式一致。
根据以上结果,我们把两个大数相乘的运算转换为了求其多项式表示形式的系数的卷积,而计算卷积最简单的方式就是通过傅里叶变换将其转变成频域的乘积,这也就是通过fft解决问题的原因。
下面给出了一个简单的python脚本实现以上算法:
1 import numpy as np 2 from scipy.fftpack import fft, ifft 3 4 5 def real_part_to_int(num): 6 """将复数的实部转换成int 7 由于精度的问题,有可能会把5表示成4.99999999999999,这时需要做一下额外处理 8 :param num: 9 :returns: 10 :rtype: 11 12 """ 13 real_part = num.real 14 lower = np.floor(real_part) 15 upper = np.ceil(real_part) 16 if upper-real_part < real_part-lower: 17 return int(upper) 18 else: 19 return int(lower) 20 21 22 ## 按照个位,十位,百位的顺序排列 23 ## 对应a=1234,b=5678 24 a = [4, 3, 2, 1, 0, 0, 0, 0] 25 b = [8, 7, 6, 5, 0, 0, 0, 0] 26 27 28 output = ifft(fft(a)*fft(b)) 29 a_prod_b = [] 30 increase = 0 31 # 从低位到高位进行进位处理 32 for ii in output: 33 ii = real_part_to_int(ii) 34 ii += increase 35 print(ii) 36 if ii > 10: 37 increase = ii//10 38 a_prod_b.append(ii % 10) 39 else: 40 increase = 0 41 a_prod_b.append(ii) 42 ## 按照正常数字顺序进行打印 43 print(a_prod_b[-1::-1])