{递归试快速傅立叶变换(FFT)乘法}
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define PI 3.1415926535897932384626

typedef struct ComplexTag {
  Real R,I;
} Complex;

long FFTLengthMax;
Complex * OmegaFFT;
Complex * ArrayFFT0, * ArrayFFT1;
Complex * ComplexCoef; 
double FFTSquareWorstError;
long AllocatedMemory;

/*初始化*/
void InitializeFFT(long MaxLength)
{
  long i;
  Real Step;

  FFTLengthMax = MaxLength;
  OmegaFFT = (Complex *) malloc(MaxLength/2*sizeof(Complex));
  ArrayFFT0 = (Complex *) malloc(MaxLength*sizeof(Complex));
  ArrayFFT1 = (Complex *) malloc(MaxLength*sizeof(Complex));
  ComplexCoef = (Complex *) malloc(MaxLength*sizeof(Complex));
  Step = 2.*PI/(double) MaxLength;
  for (i=0; 2*i<MaxLength; i++) {
    OmegaFFT[i].R = cos(Step*(double)i);
    OmegaFFT[i].I = sin(Step*(double)i);
  }
  FFTSquareWorstError=0.;
  AllocatedMemory = 7*MaxLength*sizeof(Complex)/2;
}

/*递归FFT*/
void RecursiveFFT(Complex * Coef, Complex * FFT, long Length, long Step, 
		  long Sign)
{
  long i, OmegaStep;
  Complex * FFT0, * FFT1, * Omega;
  Real tmpR, tmpI;

  if (Length==2) {
    FFT[0].R = Coef[0].R + Coef[Step].R;
    FFT[0].I = Coef[0].I + Coef[Step].I;
    FFT[1].R = Coef[0].R - Coef[Step].R;
    FFT[1].I = Coef[0].I - Coef[Step].I;
    return;
  }

  FFT0 = FFT;
  RecursiveFFT(Coef     ,FFT0,Length/2,Step*2,Sign);
  FFT1 = FFT+Length/2;
  RecursiveFFT(Coef+Step,FFT1,Length/2,Step*2,Sign);

  Omega = OmegaFFT;
  OmegaStep = FFTLengthMax/Length;
  for (i=0; 2*i<Length; i++, Omega += OmegaStep) {
    /* 递归公式:
       FFT[i]          <-  FFT0[i] + Omega*FFT1[i]
       FFT[i+Length/2] <-  FFT0[i] - Omega*FFT1[i],
       Omega = exp(2*I*PI*i/Length) */
    tmpR = Omega[0].R*FFT1[i].R-Sign*Omega[0].I*FFT1[i].I;
    tmpI = Omega[0].R*FFT1[i].I+Sign*Omega[0].I*FFT1[i].R;
    FFT1[i].R = FFT0[i].R - tmpR;
    FFT1[i].I = FFT0[i].I - tmpI;
    FFT0[i].R = FFT0[i].R + tmpR;
    FFT0[i].I = FFT0[i].I + tmpI;
  }
}

/* 计算 Coef 的复数傅立叶变换 FFT */
void FFT(Real * Coef, long Length, Complex * FFT, long NFFT)
{
  long i;
  /* 实数系数传给复数系数的实部 */
  for (i=0; i<Length; i++) {
    ComplexCoef[i].R = Coef[i];
    ComplexCoef[i].I = 0.;
  }
  for (; i<NFFT; i++)
    ComplexCoef[i].R = ComplexCoef[i].I = 0.;

  RecursiveFFT(ComplexCoef,FFT,NFFT,1,1);
}

/* 计算 FFT 的反复数傅立叶变换 Coef */
void InverseFFT(Complex * FFT, long NFFT, Real * Coef, long Length)
{
  long i;
  Real invNFFT = 1./(Real) NFFT, tmp;

  RecursiveFFT(FFT, ComplexCoef, NFFT, 1, -1);
  for (i=0; i<Length; i++) {
    /* 四舍五入 ComplexCoef[i].R/NFFT */
    tmp = invNFFT*ComplexCoef[i].R;
    Coef[i] = floor(0.5+tmp);
    if ((tmp-Coef[i])*(tmp-Coef[i])>FFTSquareWorstError)
      FFTSquareWorstError = (tmp-Coef[i])*(tmp-Coef[i]);
  }
}

void Convolution(Complex * A, Complex * B, long NFFT, Complex * C)
{
  long i;
  Real tmpR, tmpI;

  for (i=0; i<NFFT; i++) {
   tmpR = A[i].R*B[i].R-A[i].I*B[i].I;
    tmpI = A[i].R*B[i].I+A[i].I*B[i].R;
    C[i].R = tmpR;
    C[i].I = tmpI;
  }
}

void MulWithFFT(Real * ACoef, long ASize,
                Real * BCoef, long BSize,
                Real * CCoef)
{
  long NFFT = 2;

  while (NFFT<ASize+BSize)
    NFFT *= 2;

  if (NFFT>FFTLengthMax) {
    printf("Error, FFT Size is too big in MulWithFFT\n");
    return;
  }
  FFT(ACoef, ASize, ArrayFFT0, NFFT);
  FFT(BCoef, BSize, ArrayFFT1, NFFT);
  Convolution(ArrayFFT0,ArrayFFT1,NFFT,ArrayFFT0);
  InverseFFT(ArrayFFT0,NFFT,CCoef, ASize+BSize-1);
}




{递推试快速傅立叶变换(FFT)乘法}
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 :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
   {计算W^i (i=0--n)}
   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;
   {计算排序用数组, 如n=8时,lst=(0,4,2,6,1,5,3,7)}
   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 :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;
 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;
{用FFT乘法}
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 :Extended;
 Begin
   n:=1;
   logn:=0;
   while n<(aSize+bSize) do Begin n:=n shl 1; Inc(logn); End;  {确定n}
   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;
   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);
   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;
   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);
   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.