快速乘法

来源：互联网发布：邯郸市教育局网络平台编辑：程序博客网时间：2024/05/17 08:18

快速乘法

先让我们看ab乘cd是怎样进行的:

* cd

adR+bd

acR²+bcR

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

可见进行了四次乘法.

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

在上面的过程中, 似乎我们觉得这个算法并没有做任何多余的事情, 因为我们必须考虑每一组,它们的乘积对最终的结果都会产生影响. 而且我们不能通过调整计算顺序从根本上降低算法时间复杂度, 至多只能使其常数因子稍小一些.

要改进乘法速度就要设法减少乘法次数!

减少的途径, 主要有两种:

1. 二分法

2. 快速傅立叶变换

二分法

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;

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

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 次). 但是, 分成更多部分, 递归的空间消耗迅猛增加. 且算式复杂.

2. 快速傅立叶变换 (FFT)

----翻译,修改自 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 变换

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

反Fourier 变换:

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

直接按定义计算显然低效.

2.2 快速傅立叶变换(Fast Fourier Transform, FFT)

快速傅立叶变换是一个在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

2.3 FFT 乘法

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

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

计算大小为 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 时最好.

///////////////////////////////////////

涉及程序:

{$N+}{D-,E+,F-,G+,I-,L-,P+,X+,Y-}{Q-,R-,S-,T-,V-}{M 65520,0,655360}Program FFTMul;Uses dos;Const MaxSize=1 shl 4; pi=3.141592653589793238; Bas=100000; LBas=5; invBas=1/Bas;type complex=record         i,r:Extended;        end;Var i, j, k, n, asize, bSize :Longint;    WorstEr{, MaxNum }:Extended;    aCoef, bCoef, ccoef :Array[0..Pred(MaxSize)] of LongInt;    Time1,Time2,th2,tm2,ts2,tms2,th1,tm1,ts1,tms1:Word;Procedure Initial(Var w:Array of Complex; Var lst:Array of word; n:LongInt); Var i, j, n2, m :LongInt;     tmp :Extended; Begin   i:=0;   tmp:=2*pi/n;   n2:=n shr 2;   For i:=0 to n2 do Begin     w[i].i:=sin(tmp*i);     w[i].r:=Sqrt(1-sqr(w[i].i));   End;   i:=Pred(n2);   j:=succ(n2);   While i>=0 do Begin     w[j].i:= w[i].i;     w[j].r:=-w[i].r;     Dec(i);  Inc(j);   End;   n2:=n;  lst[0]:=0;  m:=1;   i:=1;   While n2<>1 do Begin     n2:=n2 shr 1;     m :=m  shl 1;     j:=0;     While i<m do Begin       lst[i]:=lst[j]+n2;       inc(i);  inc(j);     End;   End; End;Procedure PZero(x:LongInt); Begin   Inc(x);   While x<(Bas*0.1) do Begin     Write('0');     x:=x*10;    End; End;Procedure Print(a:Array of LongInt; n:LongInt); Var i :LongInt; Begin   Write (a[0]);   For i:=1 to pred(n) do begin     PZero (a[i]);     Write (a[i]);   end; End;Procedure FFT(Const a:Array of Complex; Var b:Array of Complex;Const w:Array of Complex;Const lst:Array of Word; n:Longint); Var per, m, n2, n1, j, k, i2, p :LongInt;     tmp1, tmp2 :Extended;     Tmpb, Tmpb2 :Complex; Begin   For j:=0 to pred(n) do  b[j]:=a[lst[j]];   n2:=1;  per:=1;   n1:=n shr 1;   While per<n do Begin     m:=n2;  j:=0;  k:=n2;     per:=per shl 1;     While k<n do Begin       i2:=0;       While j<m do Begin         Tmpb:=b[j];         Tmp1:=b[k].r*w[i2].r-b[k].i*w[i2].i;         Tmp2:=b[k].r*w[i2].i+b[k].i*w[i2].r;         b[j].r:=Tmpb.r+Tmp1;         b[j].i:=Tmpb.i+Tmp2;         b[k].r:=Tmpb.r-Tmp1;         b[k].i:=Tmpb.i-Tmp2;         Inc(j);   Inc(k);   Inc(i2,n1);       End;       Inc(m, per);  Inc(j, n2);  Inc(k, n2);     End;     n2:=per;     n1:=n1 shr 1;   End; End;Procedure invFFT(Const a:Array of Complex; Var b:Array of Complex;Const w:Array of Complex;Const lst:Array of Word; n:Longint); Var per, m, n1, n2, j, k, i2 :LongInt;     Tmp1, Tmp2 :Extended;     Tmpb :Complex; Begin   For j:=0 to pred(n) do  b[j]:=a[lst[j]];   n2:=1;  per:=1;   n1:=n shr 1;   While per<n do Begin     m:=n2;  j:=0;  k:=n2;     per:=per shl 1;     While k<n do Begin       i2:=0;       While j<m do Begin         Tmpb:=b[j];         Tmp1:=b[k].r*w[i2].r+b[k].i*w[i2].i;         Tmp2:=b[k].r*w[i2].i-b[k].i*w[i2].r;         b[j].r:=Tmpb.r+Tmp1;         b[j].i:=Tmpb.i-Tmp2;         b[k].r:=Tmpb.r-Tmp1;         b[k].i:=Tmpb.i+Tmp2;         Inc(j);  Inc(k);  Inc(i2,n1);       End;       Inc(m, per);  Inc(j, n2);  Inc(k, n2);     End;     n2:=per;     n1:=n1 shr 1;   End;   Tmp1:=1/n;   For j:=0 to pred(n) do     b[j].r:=b[j].r*Tmp1; End;Procedure MulFFT(Const ACoef:Array of LongInt; ASize:LongInt;Const BCoef:Array of LongInt; BSize:LongInt;                   Var CCoef:Array of LongInt); Var a, b, tmp :Array[0..MaxSize] of complex;     lst :Array[0..pred(MaxSize)] of Word;     w   :Array[0..MaxSize shr 1] of complex;     n, i :LongInt;     logn :Byte;     TmpR, num{,maxr,maxi }:Extended; Begin   n:=1;   logn:=0;   while n<(aSize+bSize) do Begin n:=n shl 1; Inc(logn); End;   Initial(w,lst,n);   for i:=0 to ASize do Begin     a[i].r:=ACoef[i];     a[i].i:=0;   End;   for i:=asize to n do Begin     a[i].r:=0;   a[i].i:=0;   End;   for i:=0 to BSize do Begin     b[i].r:=BCoef[i];     b[i].i:=0;   End;   for i:=bsize to n do Begin     b[i].r:=0;   b[i].i:=0;   End;{   maxr:=0;   maxi:=0;}   FFT(a, b, w, lst, n);   FFT(b, tmp, w,lst, n);   For i:=0 to pred(n) do Begin     a[i].r:=tmp[i].r*b[i].r - tmp[i].i*b[i].i;     a[i].i:=tmp[i].r*b[i].i + tmp[i].i*b[i].r;   End;   invFFT(a, b, w, lst, n);{   For i:=0 to pred(n) do if Abs(b[i].r)>Maxr then maxr:=Abs(b[i].r);   For i:=0 to pred(n) do if Abs(b[i].i)>Maxi then maxi:=Abs(b[i].i);}   TmpR:=0;   For i:=pred(n) downto 1 do Begin     num:=b[pred(i)].r+TmpR;     TmpR:=Int((num+0.1)*invBas);     CCoef[i]:=Round(num-TmpR*Bas);     If ABS(num-TmpR*Bas-CCoef[i]) > WorstEr Then WorstEr:=Abs((num-TmpR*Bas)-CCoef[i]);   End;   CCoef[0]:=Round(TmpR); End;Procedure SquareFFT(Const ACoef:Array of LongInt; ASize:LongInt; Var CCoef:Array of LongInt); Var a, b :Array[0..MaxSize] of complex;     lst :Array[0..pred(MaxSize)] of Word;     w   :Array[0..MaxSize shr 1] of complex;     n, i :LongInt;     logn :Byte;     TmpR, num{,maxr,maxi }:Extended; Begin   n:=1;   logn:=0;   while n<(aSize+aSize) do Begin n:=n shl 1; Inc(logn); End;   Initial(w,lst,n);   for i:=0 to Pred(ASize) do Begin     a[i].r:=ACoef[i];     a[i].i:=0;   End;   for i:=asize to n do Begin     a[i].r:=0;  a[i].i:=0;   End;{   maxr:=0;   maxi:=0;}   FFT(a, b, w, lst, n);   For i:=0 to pred(n) do Begin     a[i].r:=b[i].r*b[i].r - b[i].i*b[i].i;     a[i].i:=b[i].r*b[i].i + b[i].i*b[i].r;   End;   invFFT(a, b, w, lst, n);{   For i:=0 to pred(n) do if Abs(b[i].r)>Maxr then maxr:=Abs(b[i].r);   For i:=0 to pred(n) do if Abs(b[i].i)>Maxi then maxi:=Abs(b[i].i);}   TmpR:=0;   For i:=Pred(n) downto 1 do Begin     num :=b[Pred(i)].r+TmpR;     TmpR:=Int((num+0.001)*invBas);     CCoef[i]:=Round(num-TmpR*Bas);     If ABS(num-TmpR*Bas-CCoef[i]) > WorstEr Then WorstEr:=Abs((num-TmpR*Bas)-CCoef[i]);   End;   CCoef[0]:=Round(TmpR); End;Begin  Assign(output,'D:/Mulout.txt');  Rewrite(output);  Randomize;  WorstEr:=0;  fillchar(ACoef,sizeof(ACoef),0);  fillchar(BCoef,sizeof(BCoef),0);  asize:=Maxsize shr 1;  BSize:=Maxsize shr 1;  For i:=0 to pred(Asize) do    ACoef[i]:={Round(Bas-1)}Trunc(random*Bas);  Print(ACoef, Asize); Writeln; writeln('*');  For i:=0 to pred(Bsize) do    BCoef[i]:={Round(Bas-1)}Trunc(random*Bas);  Print(BCoef, BSize); Writeln; Writeln('='); Gettime(th1,tm1,ts1,tms1);   SquareFFT(ACoef,ASize,CCoef);{MulFFT(ACoef, ASize, BCoef, BSize, CCoef);} Gettime(th2,tm2,ts2,tms2);  Print(CCoef, Asize+BSize); Writeln;  Writeln('Worst Error is ',WorstEr);  time1:=60*(tm2-tm1)+ts2-ts1;  If tms2<tms1 then Begin time2:=100+tms2-tms1; Dec(time1); End                    else time2:=tms2-tms1;  Writeln('Used ',time1:4,'.',time2:2,'  second');  Close(output);End.