当前位置：网站首页>A concise course of fast Fourier transform FFT

A concise course of fast Fourier transform FFT

2022-04-23 09:44:00 【Zimba_】

problem

Given length is $n$ The polynomial $f(x)=\sum_{i=0}^{n-1}a_{i}x^{i}$ And length $m$ polynomial $g(x)=\sum_{i=0}^{m-1}b_{i}x^{i}$ . Solve polynomials $h(x)=f(x)\times g(x)$ .
$(1\leq n\leq 10^{5})$

Coefficient representation and point value representation

Two representations of polynomials , Coefficient representation and point value representation .
Coefficient representation ： We store $f (x)$ Each of $x^{i}$ The coefficient of $a_{i}$ .
Point value representation ： As everyone knows , $N$ A point can determine a maximum order as $N - 1$ polynomial , So we can store any number on a polynomial $N$ A little bit , To determine this polynomial .

A property of point valued representation is , For polynomials $f$ A point value of $x,y_{1})$ , Sum polynomial $g$ A point value of $x,y_{2})$ , You can get a polynomial multiplied by two polynomials $h$ The abscissa is $x$ The point at is $(x,y_{1}\times y_{2})$ .

So for the length $n$ The polynomial $f$ And length $m$ The polynomial $g$ When doing multiplication , We can figure out $f$ Of $n + m - 1$ A point value and $g$ Of $n + m - 1$ A point value , And the value of two polynomial points is required $x$ identical .
We can calculate this through traversal $n + m - 1$ The point value of the abscissa is multiplied by the ordinate , To get the polynomial after multiplication $n + m - 1$ A point value , That is, the point valued representation of the new polynomial .
( Be careful ： Although we store $f$ It only needs $n$ A point value , But in order to calculate the length $n + m - 1$ Of $h$ , We still have to get $f$ Of (n+m-1) A point value )

therefore , If we can pass $o (n l o g n)$ The way to , Convert coefficient representation to point value representation , Then take two polynomial point values $o (n)$ Multiply , Finally through $o (n l o g n)$ The method converts the point value representation back to the coefficient representation , We can solve this problem .

Later, we will introduce how to convert the coefficient representation to the point value representation .

Discrete Fourier transform （DFT）

For functions without periods , We can think of its period as $\infty$ , Use infinite periodic functions to fit .
Empathy , For a length of $N$ Discrete polynomials $x (n) (n = 0, 1, 2 . . N - 1)$ , Can pass $N$ A periodic function to fit .（ $N$ A periodic function can only fit $N$ A little bit , It cannot be used to fit continuous polynomial functions ）
From this, we get the expansion of discrete Fourier series as follows .
$X(k)=\sum_{n=0}^{N-1}x(n)e^{-2\pi ik\frac{n}{N}}$
among $e^{-2\pi ik\frac{n}{N}}=cos(-2\pi k\frac{n}{N})+isin(-2\pi k\frac{n}{N})$ It's about $n$ The period of is $N$ Function of .（ By Euler formula $e^{ix}=cos x+isin x$ have to ）

Inverse discrete Fourier transform （IDFT）

For the transformed $X (n)$ , We can get the original function through the inverse discrete Fourier transform formula $x (n)$ . And $D F T$ The formula is similar , On the original basis $i$ Turned into $- i$ , Multiply one in front $\frac{1}{N}$ . The formula is as follows .
$x(n)=\frac{1}{N}\sum_{k=0}^{N-1}X(k)e^{2\pi ik\frac{n}{N}}$

The fast Fourier transform （FFT）

As mentioned before , The core problem we need to solve is , Calculate the of polynomials $N$ A point value .
The general idea is , We find the function after discrete Fourier transform N A point value , After multiplying the new polynomial by the point value , Then it is inversely transformed back to the coefficient representation .

Let's look at the formula first $D F T$ The formula of .
In order to simplify the formula , We put $e^{\frac{2\pi ik}{N}}$ use $W_{N}$ Instead of .( $W_{N}$ It's a story about $k$ Function of )
The formula is as follows .
$X(k)=\sum_{n=0}^{N-1}x(n)W_{N}^{n}(k)$
Then observe $W_{N}$ as follows .
$W_{N}(k)=e^{\frac{2\pi i}{N}}=cos(\frac{2\pi k}{N})+isin(\frac{2\pi k}{N})$
Find out $W_{N}$ It's a cycle of $N$ Function of .

thus , We find that the transformed function $X (k)$ By $N$ The period is $N$ The power of the function of , Multiply by a factor , Add up to make . And this $N$ The real and imaginary parts of a function are trigonometric functions （ $c o s$ and $s i n$ ）.

It's easy to find us $W_{N}(k)=W_{N}(k+N)$ , But it doesn't matter to us . The more important property comes from another property of trigonometric function —— symmetry .
We put $[0,2\pi]$ In two parts $[0,\pi)$ , and $(\pi,2\pi]$ , We can see that we calculate the left $t$ Place the value of the , You can get the right $\pi +t$ Place the value of the .
To put it bluntly $sin(t)=-sin(\pi +t)$ and $cos(t)=-cos(\pi +t)$ .

Corresponding to the above, we get a conclusion , We just have to figure out $W_{N}(t)$ , You can get $W_{N}(\frac{N}{2}+t)$ . namely $W_{N}(t)=-W_{N}(\frac{N}{2}+t)$ .

The above is just a simple property derivation , The formal Fourier transform is shown below .

In short, we should divide and conquer , This formula can be organically combined with divide and conquer .（ I don't know how he came up with it , The formula is derived in 《 Digital signal processing 》 The fourth edition of P122）
We separate the odd and even terms of enumeration .
$X(k)=\sum_{n For the even }x(n)W_{N}^{n}(k)+\sum_{n It's odd }x(n)W_{N}^{n}(k) \\=\sum_{r=0}^{\frac{N}{2}-1}x(2r)W_{N}^{2r}(k)+\sum_{r=0}^{\frac{N}{2}-1}x(2r+1)W_{N}^{2r+1}(k) \\=\sum_{r=0}^{\frac{N}{2}-1}x_{1}(r)W_{N}^{2r}(k)+W_{N}(k)\sum_{r=0}^{\frac{N}{2}-1}x_{2}(r)W_{N}^{2r}(k) \\=\sum_{r=0}^{\frac{N}{2}-1}x_{1}(r)W_{\frac{N}{2}}^{r}(k)+W_{N}(k)\sum_{r=0}^{\frac{N}{2}-1}x_{2}(r)W_{\frac{N}{2}}^{r}(k) \\=X_{1}(k)+W_{N}(k)X_{2}(k)$
among $X_{1}$ yes $x_{1}$ The Fourier transform of , $X_{2}$ yes $x_{2}$ The Fourier transform of .
There is one step $W_{N}^{2r}(k)=W_{\frac{N}{2}}^{r}$ I forgot the certificate in front , Those who want to read must be able to prove it at once .

This is the place , All the preparations have been done , It is easy to find that as long as divide and conquer is done well in the point value of two small polynomials , The big ones come out naturally . So we can divide and conquer .

The final formula is as follows .
$X(k)=X_{1}(k)+W^{k}_{N}X_{2}(k),k=0,1,…,\frac{N}{2}-1 \\X(k+\frac{N}{2})=X_{1}(k)-W^{k}_{N}X_{2}(k),k=0,1,…,\frac{N}{2}-1$

But there is one detail , Calculation $W_{N}(k)$ The complexity of $o (l o g n)$ . And at every level of divide and conquer, we need to ask $N$ Time $W$ function , A total of $l o g n$ layer , Eventually, our complexity becomes $o (n l o g n l o g n)$ .
So the solution is mentioned earlier $W_{N}(t)=-W_{N}(\frac{N}{2}+t)$ , We only need to calculate half of each floor $W$ , That is to say $\frac{N}{2}$ individual $W$ , So I'm calculating $W$ The total cost on the becomes $o(\frac{n}{2}lognlogn)$ . There are pots and pans .

Then I found another place to optimize . since $W_{N}(k)=-W_{N}(\frac{N}{2}+k)$ To reduce $W$ The number of calculations （ Mainly to reduce $s i n$ and $c o s$ The number of calculations ）, It mainly comes from trigonometric function $[0,\pi]$ You can push the show $[\pi,2\pi]$ . This only needs to be calculated on each floor $\frac{N}{2}$ Time $c o s$ and $s i n$ , So let's go down one more layer , So that each layer only needs to calculate $\frac{N}{4}$ Time $c o s$ and $s i n$ .
So in the same way , Trigonometric functions know $[0,\frac{\pi}{2}]$ You can figure out $[\frac{\pi}{2},2\pi]$ .
therefore , We started pushing $W_{N}(k+\frac{N}{4})$ .
$W_{N}(k+\frac{N}{4})=cos(\frac{2\pi k}{N}+\frac{\pi}{2})+isin(\frac{2\pi k}{N}+\frac{\pi}{2})=-sin(\frac{2\pi k}{N})+icos(\frac{2\pi k}{N})$
then , We found that , We are for a certain $k$ Calculated $s i n$ and $c o s$ function , It can be used to $k+\frac{N}{4}$ , $k+\frac{N}{2}$ , $k+\frac{3N}{4}$ .
therefore , In each layer , We just need to calculate $\frac{N}{4}$ Time $W$ .
The complexity is then optimized to $o(\frac{n}{4}lognlogn)$ . There are still no eggs .

In fact, we just need to calculate $W_{N}(1)$ , Yes $W_{N}(k+1)=W_{N}(k)\times W_{N}(1)$ , This is really useful . It is proved that the geometric meaning of this complex number can be considered , Is to divide the unit circle at the center of the circle into $N$ Parts of , Taking out the first part is equivalent to a rotation matrix , Power is the number of rotations .
therefore , We just need to start with the trigonometric function to calculate $W_{N}(1)$ , Just multiply it all the way back .

And then it's done .

About the conversion back coefficient , That is, inverse transformation , It is found that the formula is similar to this , So the practice is basically the same . Remember to divide by $N$ . I won't repeat it here .

For convenience , Generally, first pull the polynomial length to 2 The power of , Then transform , Convenient recursion to the bottom layer is 1. On the other hand , To use FFT When the board , Note that the length should be 2 The power of .
Although it is divided , But for the sake of constants , It is usually written in a non recursive version .

 Recursive version 
 Recursive version constant is very large , It is strongly recommended to find a non recursive version of the code with a small constant 
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;
const double pi = acos(-1);
const int MAXN = 3000005;
struct cp
{
    
    double s,v;
    cp(double ss=0,double vv=0){
    
        s=ss;v=vv;
    }
};
cp operator * (cp a,cp b)
{
    
    return cp(a.s*b.s-a.v*b.v,a.s*b.v+a.v*b.s);
}
cp operator + (cp a,cp b)
{
    
    return cp(a.s+b.s,a.v+b.v);
}
cp operator - (cp a,cp b)
{
    
    return cp(a.s-b.s,a.v-b.v);
}
cp operator / (cp a,int b)
{
    
    return cp(a.s/b,a.v/b);
}

// Preprocessing 2 To the power of cos and sin
double coss[MAXN],sinn[MAXN];

void fft(cp *a,int len,int opt)
{
    
    if(len==1)return ;
    static cp b[MAXN];
    int mid=len>>1,mid_2=mid>>1;
    for(int i=0;i<len;i++)b[(i>>1)|(i&1?mid:0)]=a[i];
    for(int i=0;i<len;i++)a[i]=b[i];
    fft(a,mid,opt);fft(a+mid,mid,opt);
    cp uni(coss[mid],opt*sinn[mid]),w(1,0);
    for(int k=0;k<len;k++)
    {
    
        if(k<mid)b[k]=a[k]+w*a[k+mid];
        else b[k]=a[k-mid]+w*a[k];
        w=w*uni;
    }
    for(int i=0;i<len;i++)a[i]=b[i];
}
void mul(cp *a,cp *b,int len)
{
    
    fft(a,len,1);
    fft(b,len,1);
    for(int i=0;i<len;i++)a[i]=a[i]*b[i];
    fft(a,len,-1);
    for(int i=0;i<len;i++)a[i]=a[i]/len;
}
string sa,sb,sans;
cp a[MAXN],b[MAXN];
int ans[MAXN];
int main()
{
    
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cin>>sa>>sb;
    int n=sa.size(),m=sb.size();
    int len=1;
    while(len<n+m-1)coss[len]=cos(pi/len),sinn[len]=sin(pi/len),len<<=1;
    for(int i=0;i<len;i++)
    {
    
        a[i]=i<n?cp(sa[n-i-1]-'0',0):cp(0,0);
        b[i]=i<m?cp(sb[m-i-1]-'0',0):cp(0,0);
    }
    
    mul(a,b,len);
    
    for(int i=0;i<len;i++)
    {
    
        ans[i]+=round(a[i].s);
        ans[i+1]+=ans[i]/10;
        ans[i]%=10;
    }
    int flag=0;
    for(int i=len-1;i>=0;i--)
    {
    
        if(ans[i]==0&&flag==0)continue;
        else
        {
    
            flag=1;
            sans+=(char)(ans[i]+'0');
        }
    }
    cout<<sans<<endl;
    return 0;
}