FFT algorithm in the sense of modulus

Written in the previous on the FFT algorithm personal thought more important derivation, detailed introduction can refer to the FFT algorithm study notes

V[n] is a real sequence of 2N in length, and V[k] represents the 2N point DFT of the real sequence. Defines two real sequences with a length of n g[n] and H[n] as

G[N]=V[2N], h[n]=v[2n+1], 0<=n<n

The following deduction can be made

The FFT algorithm used here can achieve O (nlogn) complexity of the discrete Fourier transform and the above-mentioned relationship is closely related to the last.

Below to get to the point--the FFT in the sense of mode

Still need to know about the symmetric nature of the DfT of complex sequences and some supplementary definitions

As a result, suppose that the prime number p for the modulus is 1e8, then we can "tear" the original real sequence x[n].

Let's say we're asking for the result of X[n] convolution y[n] t[n].

Assuming that Q is a number at the sqrt (p) level, we can save the x[n]/q to the real part of the complex sequence x1[n], x[n]%q to the imaginary part of the complex sequence x1[n]. At this time, to X1[n], Y1[n] to seek the DFT, and then by X1[k]*y1[k] to get t1[k], the entire operation process can produce the largest floating point is the n*q^2 level, generally still within the acceptable range.

Next consider the process of recovering the original t[n from the convolution result {t1[k]}.

Look at the composition of T1[k]

It will be almost over by here. Found above the last line equals to the right there are four products, we can put the above four products separately, to find IDFT can recover the results of x/y_re/im convolution, and then for different products, consider the need to multiply q^2, q^1 or q^0, to restore it can be.

Serve the template

namespaceFft_mo//need to have mod (1e8~1e9 level), Upmo (A, b) definition{    Const intfft_maxn=1<< -; ConstDB pi=3.14159265358979323846264338327950288L; structCP {db A, B; CP (Doublea_=0,Doubleb_=0) {a=a_,b=b_; } CPoperator+(ConstCP&AMP;RHS)Const        {            returnCP (a+rhs.a,b+rhs.b); } CPoperator-(ConstCP&AMP;RHS)Const        {            returnCP (a-rhs.a,b-rhs.b); } CPoperator*(ConstCP&AMP;RHS)Const        {            returnCP (a*rhs.a-b*rhs.b,a*rhs.b+b*RHS.A); } CPoperator! ()Const        {            returnCP (a,-C); }}NW[FFT_MAXN+1],F[FFT_MAXN],G[FFT_MAXN],T[FFT_MAXN];//a<->f,b<->g,t<~>c    intBITREV[FFT_MAXN]; voidFft_init ()//Initialize nw[],bitrev[]    {        intL=0; while((1&LT;&LT;L)!=FFT_MAXN) l++;  for(intI=1; i<fft_maxn;i++) bitrev[i]=bitrev[i>>1]>>1| ((i&1) << (l1));  for(intI=0; i<=fft_maxn;i++) NW[I]=CP ((db) Cosl (2*pi/fft_maxn*i), (db) Sinl (2*pi/fft_maxn*i)); }        //n is guaranteed to be a power of 2 for the entire number of times//FLAG=1:DFT | FLAG=-1:IDFT    voidDFT (CP *a,intNintflag=1)        {        intD=0; while((1&LT;&LT;D) *N!=FFT_MAXN) d++;  for(intI=0; i<n;i++)if(i< (bitrev[i]>>d) Swap (A[i],a[bitrev[i]&GT;&GT;D]);//notice!         for(intL=2; l<=n;l<<=1)        {            intDel=fft_maxn/l*flag;//determines whether or not the WN changes in the complex plane clockwise or counterclockwise, and changes the spacing             for(intI=0; i<n;i+=l) {CP*le=a+i,*ri=a+i+ (l>>1); CP*w=flag==1? NW:NW+FFT_MAXN;//determine the starting point of the WN                 for(intk=0;k< (l>>1); k++) {cp ne=*ri * *W; *ri=*le-ne,*le=*le+NE; Le++,ri++,w+=del; }            }        }        if(flag!=1) for(intI=0; i<n;i++) a[i].a/=n,a[i].b/=N; }        //convo (a,n,b,m,c) a[0..n]*b[0..m], C[0..n+m]    voidConvo (LL *a,intN,ll *b,intM,LL *c) { for(intI=0; i<=n+m;i++) c[i]=0; intn=2; while(n<=n+m) n<<=1;//N is the length after the C extension         for(intI=0; i<n;i++)//expand a[],b[], deposit f[],g[], pay attention to modulo{LL AA=i<=n?a[i]:0, Bb=i<=m? B[i]:0; AA%=mod,bb%=MoD; F[i]=CP (db (aa>> the), DB (aa&32767)); G[i]=CP (db (bb>> the), DB (bb&32767));        } DFT (F,n), DFT (G,N);  for(intI=0; i<n;i++)//Restore Imaginary Part Two "product" (product specific meaning see above)        {            intJ=i? N-i:0; T[i]= ((F[i]+!f[j]) * (!g[j]-g[i]) + (!f[j]-f[i]) * (G[i]+!g[j])) *CP (0,0.25); } DFT (T,n,-1);  for(intI=0; i<=n+m;i++) Upmo (C[i], (LL (t[i].a+0.5))%mod<< the);  for(intI=0; i<n;i++)//Restore the real Part Two "product"        {            intJ=i? N-i:0; T[i]= (!f[j]-f[i]) * (!g[j]-g[i]) *CP (-0.25,0) +CP (0,0.25) * (F[i]+!f[j]) * (g[i]+!G[j]); } DFT (T,n,-1);  for(intI=0; i<=n+m;i++) Upmo (C[i],ll (t[i].a+0.5) + (LL (t[i].b+0.5)%mod<< -)); }}

Templates

This blog main reference: Digital signal processing-computer-based Method (fourth edition) [United States] Sanjit K. Mitra Yu Xiangyu translation

Reprint please specify the source http://www.cnblogs.com/Just--Do--It/p/7892254.html

Thank you for reading

FFT algorithm in the sense of modulus