快速乘法
先让我们看ab乘cd是怎样进行的:
ab
* cd
adR+bd
acR2+bcR
acR2+(ad+bc)R+bd (R是进位制, 下同.)
可见进行了四次乘法.
更一般地, 两个 n 位数相乘, 要进行 n2 次单位乘法.
在上面的过程中, 似乎我们觉得这个算法并没有做任何多余的事情, 因为我们必须考虑每一组,它们的乘积对最终的结果都会产生影响. 而且我们不能通过调整计算顺序从根本上降低算法时间复杂度, 至多只能使其常数因子稍小一些.
要改进乘法速度就要设法减少乘法次数!
减少的途径, 主要有两种:
1. 二分法
2. 快速傅立叶变换
二分法
1. 用代数式的变换可得:
ab*cd = acR2+(ad+bc)R+bd = acR2 + [(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(nlog3)≈O(n1.59).
2. 也可以将大数分为三部分:
abc*def=(ad)R4+(ae+bd)R3+(af+be+cd)R2-(bf+ce)R-cf
=adR4 + [(a-b)(e-d)+ad+be]R3 + [(a-c)(f-d)+ad+cf+be]R2 + [(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 次乘法(而不是 42=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的单位复数根Wk.即
这样一来, 变换可以在 O(n logn) 的时间里完成.
2.1 Fourier 变换
设有序列 A = (a0, a1,¼, a2n-1), 对它进行 Fourier 变换, 方法如下.
反Fourier 变换:
由于 , 所以进行正反两次Fourier 变换后, 结果应除以 2n.
直接按定义计算显然低效.
2.2 快速傅立叶变换(Fast Fourier Transform, FFT)
快速傅立叶变换是一个在O(n logn)里计算序列 A 的 Fourier Transform 的方法
由于公式:
由此在O(n logn)里计算F2n(A) 的步骤:
分A为两部分.
A0 = (a0, a2,…, a2n-2), A1 = (a1,a3,…,a2n-1).
递归,计算傅立叶变换:
C = Fn(A0) = (c0,c1,...,cn-1), D = Fn(A1) = (d0,d1,...,dn-1).
合并序列, B = (b0,...,b2n-1) = F2n(A)
bk = ck + wk dk, bn+k = ck - wk dk, 0 ≤ k < n
2.3 FFT 乘法
设 n 是2的幂. 且两个大整数 X & Y 的系数小于n
在时间 O(nlog(n)), 计算 Z = XY 在时间 O(nlog(n)), 执行下列步骤 :
计算 大小为 2n 的 傅立叶变换 X* 的序列 (xj) :
X* = (x0*,x1*,...,x2n-1*)= F2n(x0,x1,...,xn-1,0,...,0)
同样, 计算 Y*的 傅立叶变换(yj) :
Y* = (y0*,y1*,...,y2n-1*) = F2n(y0,y1,...,yn-1,0,....,0)
计算 X*, Y* 对应项的乘积到 Z*
Z* = (z0*,z1*,...,z2n-1*), zi*≡xi* yi*
计算 Z* 的反傅立叶变换 Z, 并将各项除以 2n.
然后重新整理系数zi到正常格式, 得
这么做等于计算 Z = XY.
注意在X的系数后添加高次零项使n化为2的正整数次幂, Y也是.
可见, 进行了3次大小2n的FFT, 和一组系数相乘(可忽略), 因此它的复杂度是3次大小2n的FFT.
计算平方时,只需进行2次大小2n的FFT.
让我们来看一个具体实例(用100进制): 45672=？
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 { 开始变换: bk = ck + wk dk, bn+k = ck - wk dk }
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 bk = ck + w -k dk, bn+k = ck - w -k dk }
8100, 24120, 17956, 0 { 变换完毕 }
2025, 6030, 4489, 0 { 各项除以四 }
2025
6030
4489
20857489 { 加起来 }
再用计算器验算一下, 没错吧!
2.2 数字错误
note: 由于 wk 不可能精确储存(涉及无理数), 所以会产生误差, 当误差 <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.
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
