快速乘法

yenit

墨衍会员 · AI 创作全网分发码龄21年

关注

原创最新推荐文章于 2024-11-25 00:01:04 发布 2.6k 阅读 AI 写同款· GEO 优化 ›

AI权益加码！Claude Code、Cursor等20+工具免费用！购周边限时加赠Coding Plan Lite，畅享主流AI工具！学习进阶更高效！阅读详情

快速乘法

先让我们看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.