FFT

A fast Fourier transform (FFT) is an algorithm that computes the discrete Fourier transform (DFT) of a sequence, or its inverse (IDFT). Fourier analysis converts a signal from its original domain (often time or space) to a representation in the frequency domain and vice versa. The DFT is obtained by decomposing a sequence of values into components of different frequencies.[1] This operation is useful in many fields, but computing it directly from the definition is often too slow to be practical. An FFT rapidly computes such transformations by factorizing the DFT matrix into a product of sparse (mostly zero) factors.[2] As a result, it manages to reduce the complexity of computing the DFT from O ( n 2 ) {displaystyle O(n^{2})} , which arises if one simply applies the definition of DFT, to O ( n log ⁡ n ) {displaystyle O(nlog n)} , where n {displaystyle n} is the data size. The difference in speed can be enormous, especially for long data sets where N may be in the thousands or millions. In the presence of round-off error, many FFT algorithms are much more accurate than evaluating the DFT definition directly. There are many different FFT algorithms based on a wide range of published theories, from simple complex-number arithmetic to group theory and number theory.

Fast Fourier transforms are widely used for applications in engineering, science, and mathematics. The basic ideas were popularized in 1965, but some algorithms had been derived as early as 1805.[1] In 1994, Gilbert Strang described the FFT as "the most important numerical algorithm of our lifetime",[3][4] and it was included in Top 10 Algorithms of 20th Century by the IEEE journal Computing in Science & Engineering.[5] In practice, the computation time can be reduced by several orders of magnitude in such cases, and the improvement is roughly proportional to N / {displaystyle /} log N. This huge improvement made the calculation of the DFT practical; FFTs are of great importance to a wide variety of applications, from digital signal processing and solving partial differential equations to algorithms for quick multiplication of large integers.

The best-known FFT algorithms depend upon the factorization of N, but there are FFTs with O(N log N) complexity for all N, even for prime N. Many FFT algorithms only depend on the fact that e − 2 π i / N {displaystyle e^{-2pi i/N}} is an N-th primitive root of unity, and thus can be applied to analogous transforms over any finite field, such as number-theoretic transforms. Since the inverse DFT is the same as the DFT, but with the opposite sign in the exponent and a 1/N factor, any FFT algorithm can easily be adapted for it.
FFT is to solve the following problem.

[P(x)=sum_{i=0}^{n}a_ix_i=a_0+a_1x+a_2x_2+⋯+a_nx_n ]

Complex number

A=a+ib,(a,b∈ℝ,A∈ℂ)
$ i=sqrt{-1} ( ) A_1+A_2=(a_1+i* b_1)+(a_2+i* b_2)=(a_1+a_2)+i* (b_1+b_2) $$ A_1−A_2=(a_1+i* b_1)−(a_2+i* b_2)=(a_1−a_2)+i* (b_1−b_2) $

$ A_1×A_2=(a_1+i* b_1)(a_2+i* b_2)=a_1a_2+i* 2b_1b_2+i* a_1b_2+i* a_2b_1=(a_1a_2−b_1b_2)+i* (a_1b_2+a_2b_1) $

单位根

记 $ ω_n = cos(2πn)+i* sin(2πn) $,则 n 次单位根就是

$ ω_0n,ω_1n,⋯,ω{n−1}n $ 。 通项公式, $ ω^i_n=cos(2iπn)+i* sin(2iπn) i∈[1,n] $

it has the following fact

$ ω{dk}_{dn}=ωk_n $

$ ω{k+1}_n=ωk_n* ω_n $

$ ω^{frac{n}{2}}_n=-1 $

$ ω{k+frac{n}{2}}_n=-ωk_n $
DFT

[F(x)=sum_{i=0}^{n-1}a_i* x^i ]

[F_0(x)=sum_{i=0}^{frac{n}{2}-1}a_{2i}* x^i ]

[F_1(x)=sum_{i=0}^{frac{n}{2}-1}a_{2i+1}* x^i ]

[F(x)=F_0(x)+x* F_1(x) ]

[F(ω^k_n)=F_0(ω^{2* k}n)+ω^k_n* F_1(ω^{2* k} n)=F_0(ω^k_{frac{n}{2}})+ω^k_n* F_1(ω^k_{frac{n}{2}}) ]

[F(ω^{k+frac{n}{2}}n)=(-ω^k_n)=F_0(ω^{2* k}n)-ω^k_n* F_1(ω^{2* k} n)=F_0(ω^k{frac{n}{2}})-ω^k_n* F_1(ω^k_{frac{n}{2}}) ]

分治!!!
IDFT

so as we make our $ ω^k_n $ to $ -ω^k_n $,Ans to a,we can solve it.
Better

first we consider a (n) is the smallest number that n>len1+len2 and lowerbit(n)=n 我们发现分治到边界之后的下标是原下标的二进制位翻转 so

f[0]=0;
for(int i=1;i<=n;i++)
if(i&1) f[i]=f[i-1]|(n>>1);
else f[i]=f[i>>1]>>1;

as we calculate FFT we can use the following steps

$ tmp=ω^k_n* a_{frac{n}{2}+k} $

$ a_{frac{n}{2}+k}=a_{k}-tmp ( ) a_{k}=a_{k}+tmp $

//zzd's code not writer's
#include <bits/stdc++.h>
using namespace std;
const int N=1<<20;
const double PI=acos(-1.0);
struct C{
    double r,i;
    C(){}
    C(double a,double b){r=a,i=b;}
    C operator + (C x){return C(r+x.r,i+x.i);}
    C operator - (C x){return C(r-x.r,i-x.i);}
    C operator * (C x){return C(rx.r-ix.i,rx.i+ix.r);}
}a[N],b[N],w[N];
int A,B,n,L,R[N];
void FFT(C a[],int n){
    for (int i=0;i<n;i++)
	if (R[i]>i)
	    swap(a[R[i]],a[i]);
    for (int t=n>>1,d=1;d<n;d<<=1,t>>=1)
	for (int i=0;i<n;i+=(d<<1))
	    for (int j=0;j<d;j++){
		C tmp=w[tj]a[i+j+d];
		a[i+j+d]=a[i+j]-tmp;
		a[i+j]=a[i+j]+tmp;
	    }
}
int main(){
    scanf("%d",&A);A++;
    scanf("%d",&B);B++;
    for (int i=0;i<A;i++)
	scanf("%lf",&a[i].r);
    for (int i=0;i<B;i++)
	scanf("%lf",&b[i].r);
    for (n=1,L=0;n<=A+B;n<<=1,L++);
    for (int i=0;i<n;i++){
	R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
	w[i]=C(cos(2.0iPI/n),sin(2.0iPI/n));
    }
    FFT(a,n),FFT(b,n);
    for (int i=0;i<n;i++)
	a[i]=a[i]*b[i],w[i].i=-w[i].i;
    FFT(a,n);
    A--,B--;
    for (int i=0;i<=A+B;i++)
	printf("%d ",int(a[i].r/n+0.5));
    return 0;
}

//writer's
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<iostream>
using namespace std;
const int N=1<<20;
const double pai=acos(-1.0);
int pos[N],n,m,li,L;
struct C{
    double x,y;
    C(){}
    C(double a,double b){x=a,y=b;}
    C operator + (C xx){ return C(x+xx.x,y+xx.y);}
    C operator - (C xx){ return C(x-xx.x,y-xx.y);}
    C operator * (C xx){ return C(xxx.x-yxx.y,xxx.y+yxx.x);}
}a[N],b[N],c[N],w[N];
void FFT(C a[],int n){
    for(int i=0;i<n;i++)
	if(pos[i]>i)
	    swap(a[pos[i]],a[i]);
    for(int t=n>>1,d=1;d<n;d<<=1,t>>=1)
	for(int i=0;i<n;i+=(d<<1))
	    for(int j=0;j<d;j++){
		C tmp=w[tj]a[i+j+d];
		a[i+j+d]=a[i+j]-tmp;
		a[i+j]=a[i+j]+tmp;
	    }
}
int main(){
    scanf("%d%d",&n,&m);
    n++;m++;
    for(int i=0;i<n;i++) scanf("%lf",&a[i].x);
    for(int i=0;i<m;i++) scanf("%lf",&b[i].x);
    for(li=1,L=0;li<=n+m;li<<=1,L++);
    for(int i=0;i<li;i++){
	if(i&1) pos[i]=pos[i-1]|(li>>1);
	else pos[i]=pos[i>>1]>>1;
	w[i]=C(cos(2.0ipai/li),sin(2.0ipai/li));
    }
    //for(int i=0;i<li;i++) printf("%d ",pos[i]);
    FFT(a,li);
    //for(int i=0;i<li;i++) printf("%lf %lf
",a[i].x,a[i].y);
    FFT(b,li);
    for(int i=0;i<li;i++){
	c[i]=a[i]*b[i];
	w[i].y=-w[i].y;
    }
    FFT(c,li);
    for(int i=0;i<=n+m-2;i++)
	printf("%d ",int(c[i].x/li+0.5));
    return 0;
}
原文地址:https://www.cnblogs.com/ac_forever/p/10390099.html