大数相乘是每个编程语言绕不开的话题,因为一个变量的存储空间有限,因此无法表现一个巨大的数字。今天我们要讲的就是js实现的大数相乘,先来看看结果:

快速傅立叶变换之大数相乘 gif 快速傅立叶变换之大数相乘

注:兼容小数和负数。

思路分析

最简单的大数相乘是模拟笔算,具体如下:

background Layer 1 7 3 9 5 8 * 1.新建一个数组,放入739 9 3 7 2.里面每个数字*8 72 24 56 3.进位 2 1 9 5 4.739的每个数字*5加入对应的下标 4.739的每个数字*5加入对应的下标 2 1+5*9 9+5*3 5+5*7 2 46 24 40 5.进位 2 6 8 2 4

这是最简单的实现方法,不过今天要介绍一个更效率的算法,叫快速傅立叶变换(FFT)。它可以将普通n^2复杂度降低到nlog2(n)。

基础知识

x^2=1有几个根?大家都知道,有两个,分别是+1和-1。那么x^4呢?从表面上看也是+1和-1这两个根,那是在实数界,如果在虚数界就不一定了。在虚数界里i^2=-1,因此,在虚数界,x^4有四个根,(1,0),(0,i),(-1,0),(0,-i)。

background Layer 1 (1,0) (0,i) (-1,0) (0,-i) 实部 虚部

其实x^n=1有n个复数根,它们均匀的平分一个圆。

background Layer 1 实部 虚部 (1,0) (-1/2,√3/2i) (-1/2,-√3/2i) x^3的三个复数根 实部 虚部 (1,0) (0,i) (0,-i) x^4的四个复数根 (-1,0) 实部 虚部 (1,0) (0,i) (0,-i) x^8的八个复数根 (-1,0) (√2/2,√2/2i) (-√2/2,√2/2i) (-√2/2,-√2/2i) (√2/2,-√2/2i)

大家可以发现,不管怎么平分,(1,0)是一定存在的,因为不管是几次方,1^n=1是肯定的。

我们设定x^n的n个根分别为W(n,0),W(n,1)....W(n,n-1),那么可以知道W(n,k)=cos(2*π*k/n)+sin(2*π*k/n)*i,其中(k在[0,n-1]之间),不过W(n,k)还有另外一种表示方法,即W(n,k)=(cos(2*π/n)+sin(2*π/n))^k。

background Layer 1 a (cosa,sina*i) 实部 虚部 2a (cos2*a,sin2*a*i) 实部 虚部 ix W(n,0)=(1,0) e 假设有n个根,单位旋转角度为2*π /n,那么单位复根为_W=cos(2*π /n)+sin(2*π /n)*i =cos(x)+sin(x)*i 根据欧拉公式 A B A 可知B=A*A,意味着两个相同的复数相乘得到的复数是角度翻倍。 i*a e =cos(a)+sin(a)*i 可知,A= i*2*a e =cos(2*a)+sin(2*a)*i ,B= W(n,1)=(1,0)*_W W(n,2)=(1,0)*_W*_W W(n,k)=(1,0)*_W^k ....

FFT算法

讲了这么多,是不是懵了,和我们的乘法有毛关系?不急,不急,现在就来屡屡和乘法的关系。

任何一个数字,都可以表示为多项式的样子,例如1234,可以表示为1*x^4+2*x^3+3*x^2+4(其中x=10).

那么两个数字相乘,可以看作是多项式的相乘法:

background Layer 1 12*48 (x+2)*(4*x+8) 直接相乘 (4*x^2+16*x+16) 根据特定的x得到(x+2)的y 根据特定的x得到(4*x+8)的y (y11*y21,y12*y22.....) 两个方程的y相乘 根据这些y的乘积得到最终方程式

12*48=576,上图中的4x^2+16x+16就是结果,为什么呢?因为计算出来的方程式需要按10进位,也就是(4,16,16)=>(5,7,6)。

上图中的直接相乘就是我们的传统做法,效率比较底。另外一条比较崎岖的方法就是我们今天的主角。它的大概意思是将(x+2)和(4*x+8)这两个方程式在某些特的x上的取值相乘,再将结果还原成方程式。现在不懂什么意思没关系,但我们注意到一点就是特殊的x,什么是特殊的x,特殊的x可以做什么?

数学家发现,2^n元2^n-1次方程式通过x^(2^n)=1的2^n个复数根,可以快速求出它们对应的y。

background Layer 1 2x+3 W(2,0)=(1,0) W(2,1)=(-1,0) 正常计算 y0=2*W(2,0)+3=(5,0) y1=2*W(2,1)+3=(1,0) FFT计算 2 3 1.拆分 2 3 2.蝴蝶操作 2 3 W(2,0) -1 3+2*W(2,0)=(5,0) 3-2*W(2,0)=(1,0) y0 y1

上图描述了一个比较简单的一元二次方程求y的过程,我们可以看到,在普通求值的情况下,将x^2中的两个根W(2,0),W(2,1)代入方程式,得到两个y。同时,我们采用FFT,同样可以快速得到这些y。这两种求y方式的结果是一样的,对于一个级数特别大的方程式来说,FFT的计算次数和效率明显高于前者,这就是FFT的核心思想。

background Layer 1 2x^2+x+1 W(4,0)=(1,0) W(4,1)=(0,1) 正常计算 y0=2*W(4,0)^2+W(4,0)+1=(4,0) y1=2*W(4,1)^2+W(4,1)+1=(-1,1) FFT计算 0 2 1 1 1.拆分(每次按照奇偶下标分成两拨) 0 1 2 1 2.蝴蝶操作 0 1 W(4,0) -1 y0 y1 W(4,2)=(-1,0) W(4,3)=(0,-1) y2=2*W(4,2)^2+W(4,2)+1=(2,0) y3=2*W(4,3)^2+W(4,3)+1=(-1,-1) 补齐成4元3次方程 0x^3+2x^2+x+1 0 1 2 1 2 1 W(4,0) -1 1+(1,0)*2=(3,0) 1-(1,0)*2=(-1,0) 1+(1,0)*0=(1,0) 1-(1,0)*0=(1,0) (3,0)+(1,0)*(1,0)=(4,0) (3,0)-(1,0)*(1,0)=(2,0) (-1,0)+(0,1)*(1,0)=(-1,1) W(4,0) W(4,1) (-1,0)-(0,1)*(1,0)=(-1,1) y2 y3 自上而下 自下而上 注: 复数加法:(a+bi)+(c+di)=(a+c)+(b+d)i 复数减法:(a+bi)-(c+di)=(a-c)+(b-d)i 复数乘法:(a+bi)*(c+di)=ac+bdi+adi+bci=(ac-bd)+(ad+bc)i

对于a元b次方程式,我们首先需要将它变成2^n元2^n-1次方程。FFT的第一步就是将方程式中的所有形参按照左奇数右偶数的方式依次分解,分解到只有一个形参为止。第二部就是将分解好的形参按照蝴蝶操作依次相乘相加,例如两个形参a,b那么它们上层的结果分别为b-W(n,0)*a,b+W(n,0)*a,左减右加。蝴蝶操作的层数为n。

background Layer 1 方程式:4x^7+3x^6+x^5+2x^4+9x^3+8x^2+3x+7 形参: 4 3 1 2 9 8 3 7 4 1 9 3 3 2 8 7 4 9 1 3 3 8 2 7 4 9 1 3 3 8 2 7 形参转换前的位置: 0 1 2 3 4 5 6 7 位置二进制: 000 001 010 011 100 101 110 111 转换后对应转换前的位置: 0 4 2 6 1 5 3 7 位置二进制: 000 100 010 110 001 101 011 111 二进制反转 一个形参对应的转换后的形参位置是该形参二进制反转后的二进制所代表的十进制数字

上图描绘的是FFT的第一步,也即拆分过程,讨论的是拆分前形参的位置和拆分后位置的关系。可以发现,要想知道一个形参拆分后的位置,只要将这个形参的二进制反转,得到的十进制也就是翻转后的位置。

background Layer 1 5x^7+6x^6+2x^5+x^4+2x^3+8x^2+9x+1 拆分: 5 6 2 1 2 8 9 1 5 2 2 9 6 1 8 1 5 2 2 9 6 8 1 1 5 2 2 9 6 8 1 1 自上而下 蝴蝶操作: 5 2 2 9 6 8 1 1 W(8,0) W(8,0) W(8,0) W(8,0) 1+W(8,0)*1 1-W(8,0)*1 8+W(8,0)*6 8-W(8,0)*6 9+W(8,0)*2 9-W(8,0)*2 2+W(8,0)*5 2-W(8,0)*5 t1 t2 t3 t4 t5 t6 t7 t8 t1+W(8,0)*t3 t1-W(8,0)*t3 m1 m2 m3 m4 m5 m6 m7 m8 m1+W(8,0)*m5 W(8,0) m1-W(8,0)*m5 m2+W(8,1)*m6 W(8,1) m2-W(8,1)*m6 m3+W(8,2)*m7 W(8,2) m3-W(8,2)*m7 m4+W(8,3)*m8 W(8,3) m4-W(8,3)*m8 y0 y1 y2 y3 y4 y5 y6 y7 2^n个复数根蝴蝶操作会进行n轮 每个小蝴蝶操作套路都是一样,b-W(n,k)*a,b+W(n,k)*a 自下而上 第3轮 1组 每组8个 W(8,0) t2+W(8,2)*t4 t2-W(8,2)*t4 W(8,2) t5+W(8,0)*t7 t5-W(8,0)*t7 W(8,0) t6+W(8,2)*t8 t6-W(8,2)*t8 W(8,2) W间隔1 第2轮 2组 每组4个 W间隔2 第1轮 4组 每组2个 W间隔4 需要注意的就是每次乘的W(n,k),每组蝴蝶操作的第一次肯定是乘以W(n,0),第二次不一定是W(n,1),它们之间不一定是以1递增的,递增值=2^n/(每组蝴蝶数*2) 不是 W(8,1)

上图描绘了一个2^8方程式的FFT算法,其中详细描述了蝴蝶操作的详细过程。需要注意的就是W(8,2)的位置,那里不是W(8,1)。

至此,我们就大致了解了FFT算法,它的目的就是快速的将形参(a0,a1,...a(n-1))转变成结果(y0,y1,...yn-1)的过程。专业一点的叫法就是将时域(形参)转变成频域(值)。

IFFT

IFFT是FFT的相反操作,上面我们说到,FFT其实就是将形参变成值,那么IFFT就是将值变成形参。IFFT的运算方法和FFT差不多,只是乘的时候乘的是W(n,-k),最后需要除以n。具体的可以用矩阵方法去证明,感兴趣的可以百度看看证明方法。

background Layer 1 y值:(1,0),(5,0) IFFT (5,0) (1,0) 蝴蝶操作 W(2,0) (5,0)+W(2,0)*(1,0)=(6,0) (5,0)-W(2,0)*(1,0)=(4,0) _a0 _a1 a0=_a0/2=3 a1=_a1/2=2 结果2x+3 自下而上 (1,0) (5,0) (1,0) (5,0) 拆分 (按照奇偶顺序拆分)

上面是一个简单的IFFT计算过程,大致流程和FFT类似,只是乘的是W(n,-k),并且最后要除以n。

FFT实现大数相乘

至此,所有的准备工作已经全部准备就绪了。还记得我们之前画的图吗?就是FFT计算流程的那张。现在我们跟着流程图,整体的跑一次:

background Layer 1 12*48 (x+2)*(4*x+8) 直接相乘 (4*x^2+16*x+16) 根据特定的x得到(x+2)的y 根据特定的x得到(4*x+8)的y (y11*y21,y12*y22.....) 两个方程的y相乘 根据这些y的乘积得到最终方程式 FFT计算大数相乘法 1.我们要找到两个乘数中位数较高的那个并将位数*2,找到不小于这个数的2^n。 这里就是2^2=4 为什么要这样?因为最后要通过IFFT还原,所以我们要按照结果的标准来计算。 2.计算(x+2)的FFT 0 1 W(4,0) 2+W(4,0)*0=(2,0) 2-W(4,0)*0=(2,0) 0 0 1 2 0 1 0 2 0 1 0 2 0 2 W(4,0) 1+W(4,0)*0=(1,0) 1-W(4,0)*0=(1,0) W(4,0) (2,0)+W(4,0)*(1,0)=(3,0) (2,0)-W(4,0)*(1,0)=(1,0) W(4,1) (2,0)+W(4,1)*(1,0)=(2,1) (2,0)-W(4,1)*(1,0)=(2,-1) y0 y1 y2 y3 2.计算(4x+8)的FFT 0 4 W(4,0) 8+W(4,0)*0=(8,0) 8-W(4,0)*0=(8,0) 0 0 4 8 0 4 0 8 0 4 0 8 0 8 W(4,0) 4+W(4,0)*0=(4,0) 4-W(4,0)*0=(4,0) W(4,0) (8,0)+W(4,0)*(4,0)=(12,0) (8,0)-W(4,0)*(4,0)=(4,0) W(4,1) (8,0)+W(4,1)*(4,0)=(8,4) (8,0)-W(4,1)*(4,0)=(8,-4) y0 y1 y2 y3 3.两个方程式的值相乘 y0=(x+2)的y0*(4x+8)的y0=(36,0) y1=(x+2)的y1*(4x+8)的y1=(12,16) y2=(x+2)的y2*(4x+8)的y2=(4,0) y3=(x+2)的y3*(4x+8)的y3=(12,-16) 4.IFFT (12,-16) (12,16) W(4,0) (36,0)+W(4,0)*(4,0)=(40,0) (36,0)-W(4,0)*(4,0)=(32,0) (4,0) (36,0) W(4,0) (12,16)+W(4,0)*(12,-16)=(24,0) (12,16)-W(4,0)*(12,-16)=(0,-32) W(4,0) (40,0)+W(4,0)*(24,0)=(64,0) (40,0)-W(4,0)*(24,0)=(16,0) W(4,-1) (32,0)+W(4,-1)*(0,-32)=(64,0) (32,0)-W(4,-1)*(0,32)=(0,0) _a0 _a1 _a2 _a3 (36,0) (12,16) (4,0) (12,-16) (12,16) (12,-16) (36,0) (4,0) (12,-16) (12,16) (4,0) (36,0) a0=_a0/4=(16,0) a1=_a1/4=(16,0) a2=_a2/4=(4,0) a3=_a3/4=(0,0) 方程式:4x^2+16x+16 进位得到576

补充

大家可能又能比较疑惑W(n,-k)是怎么计算的,这里给大家演示一下:

background Layer 1 实部 虚部 A B C 45 45 A点是W(8,0),B点是W(8,1),那么C 点是W(8,-1) 正k 负k 从A点逆时针走是正k,顺时针走是负k。

DEMO下载

点击下载 [0积分]一共下载0

其他文章

0
我要评论

评论

返回
×

我要评论

回复:

昵称:(昵称不超过20个字)

图片:

提交
还可以输入500个字