当前位置:网站首页>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 n n The polynomial f ( x ) = ∑ i = 0 n − 1 a i x i f(x)=\sum_{i=0}^{n-1}a_{i}x^{i} f(x)=∑i=0n−1aixi And length m m m polynomial g ( x ) = ∑ i = 0 m − 1 b i x i g(x)=\sum_{i=0}^{m-1}b_{i}x^{i} g(x)=∑i=0m−1bixi. Solve polynomials h ( x ) = f ( x ) × g ( x ) h(x)=f(x)\times g(x) h(x)=f(x)×g(x).
( 1 ≤ n ≤ 1 0 5 ) (1\leq n\leq 10^{5}) (1≤n≤105)
Coefficient representation and point value representation
Two representations of polynomials , Coefficient representation and point value representation .
Coefficient representation : We store f ( x ) f(x) f(x) Each of x i x^{i} xi The coefficient of a i a_{i} ai.
Point value representation : As everyone knows , N N N A point can determine a maximum order as N − 1 N-1 N−1 polynomial , So we can store any number on a polynomial N N N A little bit , To determine this polynomial .
A property of point valued representation is , For polynomials f f f A point value of ( x , y 1 ) (x,y_{1}) (x,y1), Sum polynomial g g g A point value of ( x , y 2 ) (x,y_{2}) (x,y2), You can get a polynomial multiplied by two polynomials h h h The abscissa is x x x The point at is ( x , y 1 × y 2 ) (x,y_{1}\times y_{2}) (x,y1×y2).
So for the length n n n The polynomial f f f And length m m m The polynomial g g g When doing multiplication , We can figure out f f f Of n + m − 1 n+m-1 n+m−1 A point value and g g g Of n + m − 1 n+m-1 n+m−1 A point value , And the value of two polynomial points is required x x x identical .
We can calculate this through traversal n + m − 1 n+m-1 n+m−1 The point value of the abscissa is multiplied by the ordinate , To get the polynomial after multiplication n + m − 1 n+m-1 n+m−1 A point value , That is, the point valued representation of the new polynomial .
( Be careful : Although we store f f f It only needs n n n A point value , But in order to calculate the length n + m − 1 n+m-1 n+m−1 Of h h h, We still have to get f f f Of (n+m-1) A point value )
therefore , If we can pass o ( n l o g n ) o(nlogn) o(nlogn) The way to , Convert coefficient representation to point value representation , Then take two polynomial point values o ( n ) o(n) o(n) Multiply , Finally through o ( n l o g n ) o(nlogn) o(nlogn) 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 N N Discrete polynomials x ( n ) ( n = 0 , 1 , 2.. N − 1 ) x(n) (n=0,1,2..N-1) x(n)(n=0,1,2..N−1), Can pass N N N A periodic function to fit .( N N N A periodic function can only fit N N 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 ) = ∑ n = 0 N − 1 x ( n ) e − 2 π i k n N X(k)=\sum_{n=0}^{N-1}x(n)e^{-2\pi ik\frac{n}{N}} X(k)=n=0∑N−1x(n)e−2πikNn
among e − 2 π i k n N = c o s ( − 2 π k n N ) + i s i n ( − 2 π k n N ) e^{-2\pi ik\frac{n}{N}}=cos(-2\pi k\frac{n}{N})+isin(-2\pi k\frac{n}{N}) e−2πikNn=cos(−2πkNn)+isin(−2πkNn) It's about n n n The period of is N N N Function of .( By Euler formula e i x = c o s x + i s i n x e^{ix}=cos x+isin x eix=cosx+isinx have to )
Inverse discrete Fourier transform (IDFT)
For the transformed X ( n ) X(n) X(n), We can get the original function through the inverse discrete Fourier transform formula x ( n ) x(n) x(n). And D F T DFT DFT The formula is similar , On the original basis i i i Turned into − i -i −i, Multiply one in front 1 N \frac{1}{N} N1. The formula is as follows .
x ( n ) = 1 N ∑ k = 0 N − 1 X ( k ) e 2 π i k n N x(n)=\frac{1}{N}\sum_{k=0}^{N-1}X(k)e^{2\pi ik\frac{n}{N}} x(n)=N1k=0∑N−1X(k)e2πikNn
The fast Fourier transform (FFT)
As mentioned before , The core problem we need to solve is , Calculate the of polynomials N N 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 DFT DFT The formula of .
In order to simplify the formula , We put e 2 π i k N e^{\frac{2\pi ik}{N}} eN2πik use W N W_{N} WN Instead of .( W N W_{N} WN It's a story about k k k Function of )
The formula is as follows .
X ( k ) = ∑ n = 0 N − 1 x ( n ) W N n ( k ) X(k)=\sum_{n=0}^{N-1}x(n)W_{N}^{n}(k) X(k)=n=0∑N−1x(n)WNn(k)
Then observe W N W_{N} WN as follows .
W N ( k ) = e 2 π i N = c o s ( 2 π k N ) + i s i n ( 2 π k N ) W_{N}(k)=e^{\frac{2\pi i}{N}}=cos(\frac{2\pi k}{N})+isin(\frac{2\pi k}{N}) WN(k)=eN2πi=cos(N2πk)+isin(N2πk)
Find out W N W_{N} WN It's a cycle of N N N Function of .
thus , We find that the transformed function X ( k ) X(k) X(k) By N N N The period is N N N The power of the function of , Multiply by a factor , Add up to make . And this N N N The real and imaginary parts of a function are trigonometric functions ( c o s cos cos and s i n sin sin).
It's easy to find us W N ( k ) = W N ( k + N ) W_{N}(k)=W_{N}(k+N) WN(k)=WN(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 π ] [0,2\pi] [0,2π] In two parts [ 0 , π ) [0,\pi) [0,π), and ( π , 2 π ] (\pi,2\pi] (π,2π], We can see that we calculate the left t t t Place the value of the , You can get the right π + t \pi +t π+t Place the value of the .
To put it bluntly s i n ( t ) = − s i n ( π + t ) sin(t)=-sin(\pi +t) sin(t)=−sin(π+t) and c o s ( t ) = − c o s ( π + t ) cos(t)=-cos(\pi +t) cos(t)=−cos(π+t).
Corresponding to the above, we get a conclusion , We just have to figure out W N ( t ) W_{N}(t) WN(t), You can get W N ( N 2 + t ) W_{N}(\frac{N}{2}+t) WN(2N+t). namely W N ( t ) = − W N ( N 2 + t ) W_{N}(t)=-W_{N}(\frac{N}{2}+t) WN(t)=−WN(2N+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 ) = ∑ n by accidentally Count x ( n ) W N n ( k ) + ∑ n by p. Count x ( n ) W N n ( k ) = ∑ r = 0 N 2 − 1 x ( 2 r ) W N 2 r ( k ) + ∑ r = 0 N 2 − 1 x ( 2 r + 1 ) W N 2 r + 1 ( k ) = ∑ r = 0 N 2 − 1 x 1 ( r ) W N 2 r ( k ) + W N ( k ) ∑ r = 0 N 2 − 1 x 2 ( r ) W N 2 r ( k ) = ∑ r = 0 N 2 − 1 x 1 ( r ) W N 2 r ( k ) + W N ( k ) ∑ r = 0 N 2 − 1 x 2 ( r ) W N 2 r ( k ) = X 1 ( k ) + W N ( k ) X 2 ( k ) 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) X(k)=n by accidentally Count ∑x(n)WNn(k)+n by p. Count ∑x(n)WNn(k)=r=0∑2N−1x(2r)WN2r(k)+r=0∑2N−1x(2r+1)WN2r+1(k)=r=0∑2N−1x1(r)WN2r(k)+WN(k)r=0∑2N−1x2(r)WN2r(k)=r=0∑2N−1x1(r)W2Nr(k)+WN(k)r=0∑2N−1x2(r)W2Nr(k)=X1(k)+WN(k)X2(k)
among X 1 X_{1} X1 yes x 1 x_{1} x1 The Fourier transform of , X 2 X_{2} X2 yes x 2 x_{2} x2 The Fourier transform of .
There is one step W N 2 r ( k ) = W N 2 r W_{N}^{2r}(k)=W_{\frac{N}{2}}^{r} WN2r(k)=W2Nr 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 N k X 2 ( k ) , k = 0 , 1 , … , N 2 − 1 X ( k + N 2 ) = X 1 ( k ) − W N k X 2 ( k ) , k = 0 , 1 , … , N 2 − 1 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 X(k)=X1(k)+WNkX2(k),k=0,1,…,2N−1X(k+2N)=X1(k)−WNkX2(k),k=0,1,…,2N−1
But there is one detail , Calculation W N ( k ) W_{N}(k) WN(k) The complexity of o ( l o g n ) o(logn) o(logn). And at every level of divide and conquer, we need to ask N N N Time W W W function , A total of l o g n logn logn layer , Eventually, our complexity becomes o ( n l o g n l o g n ) o(nlognlogn) o(nlognlogn).
So the solution is mentioned earlier W N ( t ) = − W N ( N 2 + t ) W_{N}(t)=-W_{N}(\frac{N}{2}+t) WN(t)=−WN(2N+t), We only need to calculate half of each floor W W W, That is to say N 2 \frac{N}{2} 2N individual W W W, So I'm calculating W W W The total cost on the becomes o ( n 2 l o g n l o g n ) o(\frac{n}{2}lognlogn) o(2nlognlogn). There are pots and pans .
Then I found another place to optimize . since W N ( k ) = − W N ( N 2 + k ) W_{N}(k)=-W_{N}(\frac{N}{2}+k) WN(k)=−WN(2N+k) To reduce W W W The number of calculations ( Mainly to reduce s i n sin sin and c o s cos cos The number of calculations ), It mainly comes from trigonometric function [ 0 , π ] [0,\pi] [0,π] You can push the show [ π , 2 π ] [\pi,2\pi] [π,2π]. This only needs to be calculated on each floor N 2 \frac{N}{2} 2N Time c o s cos cos and s i n sin sin, So let's go down one more layer , So that each layer only needs to calculate N 4 \frac{N}{4} 4N Time c o s cos cos and s i n sin sin.
So in the same way , Trigonometric functions know [ 0 , π 2 ] [0,\frac{\pi}{2}] [0,2π] You can figure out [ π 2 , 2 π ] [\frac{\pi}{2},2\pi] [2π,2π].
therefore , We started pushing W N ( k + N 4 ) W_{N}(k+\frac{N}{4}) WN(k+4N).
W N ( k + N 4 ) = c o s ( 2 π k N + π 2 ) + i s i n ( 2 π k N + π 2 ) = − s i n ( 2 π k N ) + i c o s ( 2 π k N ) 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}) WN(k+4N)=cos(N2πk+2π)+isin(N2πk+2π)=−sin(N2πk)+icos(N2πk)
then , We found that , We are for a certain k k k Calculated s i n sin sin and c o s cos cos function , It can be used to k + N 4 k+\frac{N}{4} k+4N, k + N 2 k+\frac{N}{2} k+2N, k + 3 N 4 k+\frac{3N}{4} k+43N.
therefore , In each layer , We just need to calculate N 4 \frac{N}{4} 4N Time W W W.
The complexity is then optimized to o ( n 4 l o g n l o g n ) o(\frac{n}{4}lognlogn) o(4nlognlogn). There are still no eggs .
In fact, we just need to calculate W N ( 1 ) W_{N}(1) WN(1), Yes W N ( k + 1 ) = W N ( k ) × W N ( 1 ) W_{N}(k+1)=W_{N}(k)\times W_{N}(1) WN(k+1)=WN(k)×WN(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 N 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 ) W_{N}(1) WN(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 N 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;
}
Non recursive version
版权声明
本文为[Zimba_]所创,转载请带上原文链接,感谢
https://yzsam.com/2022/04/202204230621424563.html
边栏推荐
- Flutter's loading animation is more interesting
- Exclusive thoughts and cases of JS
- Simple understanding of arguments in JS
- Three ways to create objects in JS
- Learn FPGA (from Verilog to HLS)
- Amazon cloud technology entry Resource Center, easy access to the cloud from 0 to 1
- 成功的DevOps Leader 应该清楚的3个挑战
- 1D / 1D dynamic programming learning summary
- Unfortunately, I broke the leader's confidential documents and spit blood to share the code skills of backup files
- Chinese Remainder Theorem and extended Chinese remainder theorem that can be understood by Aunt Baojie
猜你喜欢
Go language learning notes - array | go language from scratch
Leetcode题库78. 子集(递归 c实现)
Kernel PWN learning (3) -- ret2user & kernel ROP & qwb2018 core
DVWA range practice
AQS & reentrantlock implementation principle
Using JS to realize a thousandth bit
Flink 流批一体在小米的实践
kettle庖丁解牛第14篇之JSON输入
SAP pi / PO soap2proxy consumption external WS example
Leetcode0587. 安装栅栏(difficult)
随机推荐
(Extended) bsgs and higher order congruence equation
1D / 1D dynamic programming learning summary
kernel-pwn学习(3)--ret2user&&kernel ROP&&QWB2018-core
653. Sum of two IV - input BST
LeetCode 1611. The minimum number of operations to make an integer 0
JSON input of Chapter 14 of kettle paoding jieniu
Integral function and Dirichlet convolution
MySQL of database -- Fundamentals
Acquisition of DOM learning elements JS
Yyds dry goods inventory ubuntu18 0.4 install MySQL and solve error 1698: access denied for user ''root' '@' 'localhost' '
Simple understanding of arguments in JS
论文阅读《Integrity Monitoring Techniques for Vision Navigation Systems》
How to obtain geographical location based on photos and how to prevent photos from leaking geographical location
【读书笔记】《Verilog数字系统设计教程》 第5章 条件语句、循环语句和块语句(附思考题答案)
Compile and debug mysql8 with clion under MacOS x
[ACM-ICPC 2018 Shenyang Network preliminaries] J. Ka Chang (block + DFS sequence)
JS and how to judge custom attributes in H5
Using sqlmap injection to obtain the account and password of the website administrator
Easy to understand subset DP
Leetcode0587. Install fence