返回主页<-

此页的.DOC版本

涉及程序-> FFTmul

ab

*    cd  

adR+bd

acR²+bcR             

acR²+(ad+bc)R+bd        (R是进位制, 下同.)

可见进行了四次乘法.

更一般地, 两个 n 位数相乘, 要进行 n² 次单位乘法.

1.1 用代数式的变换可得:

   ab*cd = acR²+(ad+bc)R+bd = acR² + [(a-d)(b-c)+ac+bd]R + bd

可见只进行了3次乘法.

由此联想到用递归来做, 得下示算法 (n 应是2的正整数次幂):

function MULT(X, Y, n); { X 和 Y 为 2 个小于 2n 的整数, 返回结果为 X 和 Y 的乘积 XY }
begin
  if n=1 then
     if (X=1)and(Y=1) then return(S)
                      else return(0)
     else begin
        A:=X的左边n/2位;
        B:=X的右边n/2位;
        C:=Y的左边n/2位;
        D:=Y的右边n/2位;
        ml:=MULT(A,C,n/2);
        m2:=MULT(A-B,D-C,n/2);
        m3:=MULT(B,D,n/2);
        S:=S*(m1*2n+(m1+m2+m3)*2n/2+m3);
       return(S);
     end;
end;

它的时间复杂度是 O(n^log3)≈O(n^1.59).

2.2 也可以将大数分为三部分:

abc*def=(ad)R⁴+(ae+bd)R³+(af+be+cd)R²-(bf+ce)R-cf

=adR⁴ + [(a-b)(e-d)+ad+be]R³ + [(a-c)(f-d)+ad+cf+be]R² + [(b-c)(f-e)+be+cf]R + cf

  分成更多部分也行, 例如 4 部分:

abcd*efgh=aeR6+(af+be)R5+(ag+bf+ce)R4+(ah+bg+cf+de)R3+(bh+cg+df)R2+(ch+dg)R+dh

=aeR6 + [(a-b)(f-e)+ae+bf]R5+[(a-c)(g-e)+ae+cg+be]R4 +

[(a-b+c-d)(h-g+f-e)-(af+be)+(ag+bf+ce)+(bh+cg+df)-(ch+dg)+ae+dh]R3 +

[(b-d)(h-f)+bf+dh+cg]R2 + [(c-d)(h-g)+cg+dh]R + dh

观察发现斜体字部分正好是 R5, R4, R2, R 项的系数(画线部分).这对优化有很大帮助, 使得实际上只需进行 9 次乘法(而不是 4²=16 次).

但是, 分成更多部分, 递归的空间消耗迅猛增加. 且算式复杂.

----翻译,修改自 Mathematical Constants and computation  的

FFT based multiplication of large numbers

两个大整数可写成如下形式:

  X=P(B),  Y=Q(B)        B是基数 (通常 B = 10 或 10 的正整数次方) 可以看作 P , Q  是两个多项式:

有它们的乘积 R(z), 那么 XY=R(B). 这样, 问题转为求两个多项式的乘积.而一个不超过 m 次的多项式可以用 m 个点 ( 当z取 m 个不同值时的 R(z) ) 来表示. 并且适当的点值表示的多项式相乘可以在 O(n) 的时间里完成. 问题是怎样把多项式快速转为点值表示法. 

m个z应取1的单位复数根W_k.即

这样一来, 变换可以在 O(n logn) 的时间里完成.

2.1 Fourier 变换  {程序 *.pas}

 设有序列 A = (a₀, a₁,¼, a_2n-1), 对它进行 Fourier 变换, 方法如下.

 反Fourier 变换:

由于  , 所以进行正反两次Fourier 变换后, 结果应除以 2n.

直接按定义计算显然低效. {用FFT做乘法的具体步骤见 2. 用FFT做乘法 }

2.2  快速傅立叶变换(Fast Fourier Transform)   {程序 *.c} 

快速傅立叶变换是一个在O(n logn)里计算序列 A  的 Fourier Transform 的方法 

由于公式:

由此在 O(n logn)里计算F_2n(A) 的步骤 :

·              分 A 为两部分.

            A₀ = (a₀, a₂,…, a_2n-2),    &    A₁ = (a₁,a₃,…,a_2n-1).

·             递归,计算傅立叶变换:

            C = F_n(A₀) = (c₀,c₁,...,c_n-1)     &   D = F_n(A₁) = (d₀,d₁,...,d_n-1).

·             合并序列, B = (b₀,...,b_2n-1) = F_2n(A)

b_k = c_k + w^k d_k,    b_n+k = c_k - w^k d_k,    0 ≤ k < n

{ 用FFT做乘法的具体步骤见 下 }

2.3  FFT 乘法

  设 n 是2的幂. 且两个大整数 X & Y 的系数小于n

在时间 O(nlog(n)) 计算 Z = XY, 执行下列步骤 :

计算 大小为 2n 的 傅立叶变换 X^* 的序列 (x_j) :

     X^* = (x₀^*,x₁^*,...,x_2n-1^*)= F_2n(x₀,x₁,...,x_n-1,0,...,0)

同样, 计算 Y^*的 傅立叶变换(y_j) :

     Y^* = (y₀^*,y₁^*,...,y_2n-1^*) = F_2n(y₀,y₁,...,y_n-1,0,....,0)

计算 X*, Y* 对应项的乘积到 Z^*

 Z^* = (z₀^*,z₁^*,...,z_2n-1^*),      z_i^*≡x_i^* y_i^*

计算 Z* 的反傅立叶变换 Z, 并将各项除以 2n.

然后重新整理系数 z_i 到正常格式, 得

这么做等于计算 Z = XY.

注意在 X 的系数后添加高次零项使n化为2的正整数次幂, Y也是.

可见, 进行了 3 次大小 2n 的FFT,和一组系数相乘(可忽略), 因此它的复杂度是3次大小 2n 的FFT.

计算平方时,只需进行 2 次大小 2n 的FFT.

让我们来看一个具体实例 (用100进制) :  4567²=？

     0   1   2   3

A =(45  67  00  00)  { 补零 }

45,00  67,00       { 分成奇,偶两组 }

 0   1   2   3

45  00  67  00       { 每组再分成奇,偶两组.分组完成 }

    0      1

45,45  67,67         { 开始变换: b_k = c_k + w^k d_k,    b_n+k = c_k - w^k d_k }

112, 45+67i ,-22, 45-67i { 变换完毕 }

12544, -2464+6030i, 484, -2464-6030i  {各项平方}

12544, 484   -2464+6030i ,-2464-6030i {分成奇,偶两组}

12544   484   -2464+6030i    -2464-6030i  { 每组再分成奇,偶两组.分组完成 }

13028, 12060   -4928, 12060i  { 反FFT b_k = c_k + w ^-k d_k,    b_n+k = c_k - w ^-k d_k }

8100, 24120, 17956,  0    { 变换完毕 }

2025,  6030,  4489,  0      { 各项除以四 }

2025

  6030

    4489

20857489              { 加起来 }

再用计算器验算一下, 没错吧!

2.2 数字错误

note: 由于 w^k 不可能精确储存(涉及无理数), 所以会产生误差, 当误差 <0.5 时, 可以容忍, <0.25 时最好.