[Transfer] Very simple and easy to understand FFT (Fast Fourier Transform)

FFT Preface

Fast Fourier transform is the collective name for efficient and fast calculation methods that use computers to calculate the Discrete Fourier Transform (DFT), referred to as FFT. Fast Fourier transform was proposed by J.W. Cooley and T.W. Tukey in 1965. Using this algorithm can greatly reduce the number of multiplications required by the computer to calculate the discrete Fourier transform. In particular, the more sampling points N are transformed, the more significant the calculation savings of the FFT algorithm will be.

FFT (Fast Fourier Transformation) is a fast algorithm for discrete Fourier transform (DFT). That is the fast Fourier transform. It is obtained by improving the discrete Fourier transform algorithm based on the odd, even, imaginary, real and other characteristics of the discrete Fourier transform.

–Baidu Encyclopedia

FFT (Fast Fourier Transformation), Chinese name Fast Fourier Transform, used to accelerate polynomial multiplication
Naive high-precision multiplication time is O(n2)O(n2) O(n^2) O(n2), but FFTFFT FFT FFT can be O(nlog2n)O(nlog?2n) O(n\log_2 n) O(nlog2 ?n) time to solve
FFTFFT FFT FFT has a very fancy name and is difficult to understand. Other tutorials are too difficult to understand, so I just wrote it down by myself.

• It is recommended that you have some knowledge of complex numbers and trigonometric functions (it doesn’t matter if you don’t know)

• I will steal the following difficult points from the Internet

Polynomial coefficient representation and point value representation

• FFTFFT FFT FFT is actually a method that uses O(nlog2n)O(nlog?2n) O(n\log_2n) O(nlog2?n) time to convert a polynomial represented by coefficients into its Point value representation algorithm

• Polynomial coefficient representation and point value representation can beconverted

Coefficient representation

A n-1 degree n term polynomial f(x)f(x) f(x) f(x) can be expressed as f(x)=∑n?1i=0aixif(x)=∑ i=0n?1aixi f(x)=\sum^{n-1}_{i=0}a_ix^i f(x)=∑i=0n?1?ai?xi
You can also use the coefficient of each term to represent f(x)f(x) f(x) f(x), that is, f(x)={a0,a1,a2,… ,an?1}f(x)={a0,a1,a2,…,an?1} f(x)=\{a_0,a_1,a_2,…,a_{n-1} \ \} f(x)={a0?,a1?,a2?,…,an?1?}
This is the coefficient representation, which is the method usually used in mathematics classes

Point value notation

• Put the polynomial into the plane rectangular coordinate system and regard it as afunction

• Substituting nn n n different xx x x, you will get nn n n different yy y y, which are nn n n different points in the coordinate system

• Then these nn n n points uniquely determine the polynomial, that is, there is only one polynomial satisfying?k,f(xk)=yk?k,f(xk) =yk ?k,f(x_k)=y_k ?k,f(xk?)=yk?

• The reason is very simple. Connect n n n equations to form an n-element system of equations with n equations. The coefficients of each term can be solved.

Then f(x)f(x) f(x) f(x) can also be used as f(x)={(x0,f(x0)),(x1,f(x1)),(x2,f(x2 )),…,(xn?1,f(xn?1))}f(x)={(x0,f(x0)),(x1,f(x1)),(x2,f(x2 )),…,(xn?1,f(xn?1))} f(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f (x_2)),…,(x_{n-1},f(x_{n-1}))\} f(x)={(x0?,f(x0?)),(x1? ,f(x1?)),(x2?,f(x2?)),…,(xn?1?,f(xn?1?))} to represent
This isPoint Value Notation

The difference between two polynomial representations in high-precision multiplication

For two polynomials represented bycoefficients, we multiply them
Let the two polynomials be A(x),B(x)A(x),B(x) A(x),B(x) A(x),B(x)
We want to enumerate the coefficients of each bit of AA A A multiplied by the coefficients of each bit of BB B B
Then the coefficient representation does polynomial multiplicationTime complexityO(n2)O(n2) O(n^2) O(n2)

But the multiplication of two polynomials represented by point values only takes O(n)O(n) O(n) O(n) time.

What does it mean?

Let the two point-valued polynomials be f(x)={(x0,f(x0)),(x1,f(x1)),(x2,f(x2)),…,(xn?1, f(xn?1))}f(x)={(x0,f(x0)),(x1,f(x1)),(x2,f(x2)),…,(xn?1, f(xn?1))} f(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),…,(x_{ n-1},f(x_{n-1}))\} f(x)={(x0?,f(x0?)),(x1?,f(x1?)),(x2?, f(x2?)),…,(xn?1?,f(xn?1?))}g(x)={(x0,g(x0)),(x1,g(x1)), (x2,g(x2)),…,(xn?1,g(xn?1))}g(x)={(x0,g(x0)),(x1,g(x1)), (x2,g(x2)),…,(xn?1,g(xn?1))} g(x)=\{(x_0,g(x_0)),(x_1,g(x_1) ),(x_2,g(x_2)),…,(x_{n-1},g(x_{n-1}))\} g(x)={(x0?,g(x0? )),(x1?,g(x1?)),(x2?,g(x2?)),…,(xn?1?,g(xn?1?))}
Suppose their product is h(x)h(x) h(x) h(x), then h(x)={(x0,f(x0)?g(x0)),(x1,f(x1) ?g(x1)),…,(xn?1,f(xn?1)?g(xn?1))}h(x)={(x0,f(x0)?g(x0)) ,(x1,f(x1)?g(x1)),…,(xn?1,f(xn?1)?g(xn?1))} h(x)=\{(x_0, f(x_0)·g(x_0)),(x_1,f(x_1)·g(x_1)),…,(x_{n-1},f(x_{n-1})·g(x_ {n-1}))\} h(x)={(x0?,f(x0?)?g(x0?)),(x1?,f(x1?)?g(x1?)), …,(xn?1?,f(xn?1?)?g(xn?1?))}

So the time complexity here is only O(n)O(n) O(n) O(n) for an enumeration

• Suddenly it feels like high-precision multiplication can O(n)O(n) O(n) O(n)Extremelya lot of questions?

• However, the algorithm for converting from simple coefficient representation to point value representation is still O(n2)O(n2) O(n^2) O(n2), and the inverse operation is similar.

• The algorithm for converting naive coefficients into point values is called DFT (Discrete Fourier Transform), and the algorithm for converting point values into coefficients is called IDFT (Inverse Discrete Fourier Transform) strong>

• Could it be that high-precision multiplication can only be O(n2)O(n2) O(n^2) O(n2)?

DFT prerequisite knowledge & skills

Plural numbers

After all, it’s in high school, so I won’t say much.

We call numbers of the form a + bi (a and b are both real numbers) called complex numbers, where a is called the real part, b is called the imaginary part, and i is called the imaginary unit. When the imaginary part is equal to zero, this complex number can be regarded as a real number; when the imaginary part of z is not equal to zero and the real part is equal to zero, z is often called a pure imaginary number. The complex number field is the algebraic closure of the real number field, that is, any polynomial with complex coefficients always has roots in the complex number field. Complex numbers were first introduced by Cardan, a scholar in Milan, Italy, in the 16th century. Through the work of D’Alembert, De Moivre, Euler, Gauss and others, this concept was gradually accepted by mathematicians.
–Baidu Encyclopedia

The junior high school math teacher will tell you that there is no ?1?√?1 \sqrt{-1} ?1?, but only RR R R
Extended to the complex number set CC C C, define i2=?1i2=?1 i^2=-1 i2=?1, a complex number zz z z can be expressed as z=a + bi(a,b∈R)z=a + bi (a,b∈R) z=a + bi(a,b\in R) z=a + bi(a,b∈R)
Among them, aa a a is called the real part, bb b b is called the imaginary part, and ii i i is called the imaginary unit

• In the set of complex numbers, ii i i can be used to represent the square root of a negative number, such as ?7?√=7–√i?7=7i \sqrt{-7}=\sqrt{7}i ?7?=7? i

Complex numbers can also be viewed as a point on acomplex plane rectangular coordinate system, such as the following

xx x

Theplural number represented by this point (2,3)(2,3) (2,3) (2,3) is 2 + 3i2 + 3i 2 + 3i 2 + 3i, or imagine what it represents Vector is (2,3)(2,3) (2,3) (2,3)
In fact, we can also imagine (actually there is no such expression) It can be expressed as (13√,θ)(13,θ) (\sqrt{13},\theta) (13?,θ)
The module of a complex number zz z z is defined as its distance from the origin, recorded as ∣z∣=a2 + b2√∣z∣=a2 + b2 |z|=\ sqrt{a^2 + b^2} ∣z∣=a2 + b2?
The complex conjugate of a complex number z=a + biz=a + bi z=a + bi z=a + bi is a?bia?bi a-bi a?bi (the negation of the imaginary part) , recorded as zˉ=a?biz￣=a?bi \overline{z}=a-bi z=a?bi

Complex number operations

Complex numbers are not like points or vectors. They can perform four arithmetic operations just like real numbers.
Suppose the two complex numbers are z1=a + bi,z2=c + diz1=a + bi,z2=c + di z_1=a + bi,z_2=c + di z1?=a + bi,z2?=c + di, then
z1 + z2=(a + c) + (b + d)iz1 + z2=(a + c) + (b + d)i z_1 + z_2=(a + c) + (b + d)i z1? + z2?=(a + c) + (b + d)iz1z2=(ac?bd) + (ad + bc)iz1z2=(ac?bd) + (ad + bc)i z_1z_2=(ac?bd) + ( ad + bc)i z1?z2?=(ac?bd) + (ad + bc)i

The addition of complex numbers also satisfies the parallelogram rule

This picture was stolen from the Internet

That is, AB + AD=ACAB + AD=AC AB + AD=AC AB + AD=AC

Multiplication of complex numbers also has a noteworthy little property
(a1,θ1)?(a2,θ2)=(a1a2,θ1 + θ2)(a1,θ1)?(a2,θ2)=(a1a2,θ1 + θ2) (a_1,\theta_1)*(a_2,\ \theta_2)=(a_1a_2,\theta_1 + \theta_2) (a1?,θ1?)?(a2?,θ2?)=(a1?a2?,θ1? + θ2?)
That is, the module lengths are multiplied and the polar angles are added

DFT (Discrete Fourier Transform)

• Be sure to note that from here on allnn n ndefaults to integer powers of 22 2 2

For any coefficient polynomial conversion point value, of course you can take any nn n n xx x x values and substitute them into the calculation.
But violent calculation x0k,x1k,…,xn?1k(k∈[0,n))xk0,xk1,…,xkn?1(k∈[0,n)) x_k^0,x_k^1,…,x_k^{n-1}(k\in[0,n)) xk0?,xk1?,…,xkn?1?(k∈[0, n)) is of course O(n2)O(n2) O(n^2) O(n2) time
In fact, you can substitute a set of magical xx x x, and you don’t need to do so many power operations after substitution.
Of course these xx x x are not taken randomly, and taking these xx x x values should be Fourier’s idea

Think about it, if we substitute some xx x x so that the power of each xx x x equals 11 1 1, we don’t have to do all the power operations
±1±1 ±1 ±1 is okay, and if you consider imaginary numbers, ±i±i ±i ±i is also okay, but these four numbers are far from enough.

• Fourier said: Every point on this circle can be done

With the origin as the center, draw a unit circle with a radius of 11 1 1
Then all points on the unit circle can be raised to several powers to get 11 1 1
Fourier said that it needs to be equally divided into n n n. For example, n=8n=8 n=8 n=8

The orange point is the point to be taken when n=8n=8 n=8 n=8, starting from the (1,0)(1,0) (1,0) (1,0) point, counterclockwise from 00 0 Numbering starts from No. 0 and goes to No. 77 7 7
Note that the complex value represented by the point numbered kk k k is ωknωnk \omega_n^k ωnk?, then by multiplying the module length and adding the polar angle, we can know that (ω1n)k=ωkn(ωn1)k =ωnk (\omega_n^1)^k=\omega_n^k (ωn1?)k=ωnk?
Among them, ω1nωn1 \omega_n^1 ωn1? is called nn n nsubunit root, and each ωω \omega ω can be found (I am not good at trigonometric functions)
ωkn=coskn2π + isinkn2πωnk=cos?kn2π + isin?kn2π \omega_n^k=\cos{k\over n}2π + i\sin{k\over n} 2π ωnk?=cosnk?2π + isinnk?2π

In other words, this formula can also be explained like this

Note that sin2θ + cos2θ=1 sin2θ + cos2θ=1 sin^2\theta + cos^2\theta=1 sin2θ + cos2θ=1 or something, it is easy to understand

Then ω0n,ω1n,…,ωn?1nωn0,ωn1,…,ωnn?1 \omega^0_n,\omega^1_n,…,\omega^{n-1}_n ωn0? ,ωn1?,…,ωnn?1? is the x0,x1,…,xn?1×0,x1,…,xn?1 x_0,x_1,…,x_{n we want to substitute -1} x0?,x1?,…,xn?1?

Some properties of unit roots

In the process of pushing FFTFFT FFT FFT, some properties of ωω \omega ω need to be used

ωkn=ω2k2nωnk=ω2n2k \omega^k_n=\omega^{2k}_{2n} ωnk?=ω2n2k?

• The points (or vectors) they representare the same

• Proof

• ωkn2π \over 2n}2π + isin{2k\over 2n} 2π=\omega^{2k}_{2n} ωnk?=cosnk?2π + isinnk?2π=cos2n2k?2π + isin2n2k?2π=ω2n2k?

ωk + n2n=?ωknωnk + n2=?ωnk \omega^{k + {n \over 2}}_n=-\omega_n^k ωnk + 2n=?ωnk?

• The points they represent aresymmetric about the origin, the complex numbers they representhave opposite real parts, and the vectors they representare opposite in magnitude

• Proof

• ωn2n=cosn2n2π + isinn2n2π=cosπ + isinπ=?1ωnn2=cosn2n2π + isinn2n2π=cosπ + isinπ=?1 \omega^{n\over 2}_n=cos{{n\over 2}\over n} 2\pi + isin{{n\over 2}\over n}2\pi=cos\pi + isin\pi=-1 ωn2n=cosn2n2π + isinn2n2π= cosπ + isinπ=?1

• (This thing is the same as eix=cosx + isinxeix=cosx + isinx e^{ix}=cosx + isinx eix=cosx + isinx and eiπ + 1=0eiπ + 1=0 e^{i\pi} + 1=0 eiπ + 1=0 has something to do with it, I won’t talk about it if I don’t know it)

ω0n=ωnnωn0=ωnn \omega^0_n=\omega^n_n ωn0?=ωnn?

• are equal to 11 1 1, or 1 + 0i1 + 0i 1 + 0i 1 + 0i

FFT (Fast Fourier Transform)

Although DFTDFT DFT DFT has produced a bunch of very cool ωω \omega ω as xx x x values that are substituted into the polynomial
But just substituting this kind of special xx

• DFT is still violentO(n2)O(n2) O(n^2) O(n2)…

But DFTDFT DFT DFT can be done by divide and conquer, so FFT (Fast Fourier Transform) came out
First, let’s assume a polynomial A(x)A(x) A(x) A(x)
A(x)=∑n?1i=0aixi=a0 + a1x + a2x2 + … + an?1xn?1A(x)=∑i=0n?1aixi=a0 + a1x + a2x2 + … + an? 1xn?1 A(x)=\sum^{n-1}_{i=0}a_ix^i=a_0 + a_1x + a_2x^2 + … + a_{n-1}x^{n- 1} A(x)=i=0∑n?1?ai?xi=a0? + a1?x + a2?x2 + … + an?1?xn?1

Divide A(x)A(x) A(x) A(x) into two halves according to the parity of the subscripts of A(x)A(x) A(x) A(x). The right side Mention one more xx x x

A(x)=(a0 + a2x2 + … + an?2xn?2) + (a1x + a3x3 + … + an?1xn?1)A(x)=(a0 + a2x2 + … + an?2xn?2) + (a1x + a3x3 + … + an?1xn?1) A(x)=(a_0 + a_2x^2 + … + a_{n-2}x^{n-2 }) + (a_1x + a_3x^3 + … + a_{n-1}x^{n-1}) A(x)=(a0? + a2?x2 + … + an?2?xn ?2) + (a1?x + a3?x3 + … + an?1?xn?1)

=(a0 + a2x2 + … + an?2xn?2) + x(a1 + a3x2 + … + an?1xn?2)=(a0 + a2x2 + … + an?2xn?2) + x(a1 + a3x2 + … + an?1xn?2) =(a_0 + a_2x^2 + … + a_{n-2}x^{n-2}) + x(a_1 + a_3x^2 + … + a_{n-1}x^{n-2}) =(a0? + a2?x2 + … + an?2?xn?2) + x(a1? + a3?x2 + … + an?1?xn?2)

Both sides seem to be very similar, so let’s set up two more polynomials A1(x),A2(x)A1(x),A2(x) A_1(x),A_2(x) A1?(x),A2?(x), make

A1(x)=a0 + a2x + a4x2 + … + an?2xn2?1A1(x)=a0 + a2x + a4x2 + … + an?2xn2?1 A_1(x)=a_0 + a_2x + a_4x^ 2 + … + a_{n-2}x^{{n\over 2}-1} A1?(x)=a0? + a2?x + a4?x2 + … + an?2? x2n1A2(x)=a1 + a3x + a5x2 + … + an?1xn2?1A2(x)=a1 + a3x + a5x2 + … + an?1xn2?1 A_2(x)=a_1 + a_3x + a_5x^2 + … + a_{n-1}x^{{n \over 2}-1} A2?(x)=a1? + a3?x + a5?x2 + … + an ?1?x2n1

Comparison It is obvious that
A(x)=A1(x2) + xA2(x2)A(x)=A1(x2) + xA2(x2) A(x)=A_1(x^2) + xA_2(x^2) A(x) =A1?(x2) + xA2?(x2)

Suppose k(In the next few transformation steps, you need to think more about the properties of the unit root)

A(ωkn)=A1((ωkn)2) + ωknA2((ωkn)2)A(ωnk)=A1((ωnk)2) + ωnkA2((ωnk)2) A(\omega^k_n)=A_1 ((\omega^k_n)^2) + \omega^k_nA_2((\omega^k_n)^2) A(ωnk?)=A1?((ωnk?)2) + ωnk?A2?(( ωnk?)2)=A1(ω2kn) + ωknA2(ω2kn)=A1(ωkn2) + ωknA2(ωkn2)=A1(ωn2k) + ωnkA2(ωn2k)=A1(ωn2k) + ωnkA2(ωn2k) =A_1(\ omega^{2k}_n) + \omega^k_nA_2(\omega^{2k}_n)=A_1(\omega^k_{n\over2}) + \omega^k_nA_2(\omega^k_ {n\over 2}) =A1?(ωn2k?) + ωnk?A2?(ωn2k?)=A1?(ω2n?k?) + ωnk?A2?(ω2n?k?)

So for A(ωk + n2n)A(ωnk + n2) A(\omega^{k + {n\over2}}_n) A(ωnk + 2n), substitute ωk + n2nωnk + n2 \ omega^{k + {n \over 2}}_n ωnk + 2n
A(ωk + n2n)=A1(ω2k + nn) + ωk + n2nA2(ω2k + nn)A(ωnk + n2)=A1(ωn2k + n) + ωnk + n2A2(ωn2k + n) A(\omega^ {k + {n\over 2}}_n)=A_1(\omega^{2k + n}_n) + \omega^{k + {n\over 2}}_nA_2(\omega^{ 2k + n}_n) A(ωnk + 2n)=A1?(ωn2k + n?) + ωnk + 2nA2?(ωn2k + n?)=A1(ω2knωnn)?ωknA2(ω2knωnn)=A1( ωn2kωnn)?ωnkA2(ωn2kωnn) =A_1(\omega^{2k}_n\omega^n_n)-\omega^k_nA_2(\omega^{2k}_n\omega^n_n) =A1?(ωn2k ?ωnn?)?ωnk?A2?(ωn2k?ωnn?)=A1(ω2kn)?ωknA2(ω2kn)=A1(ωkn2)?ωknA2(ωkn2)=A1(ωn2k)?ωnkA2(ωn2k)=A1(ωn2k) ?ωnkA2(ωn2k) =A_1(\omega^{2k}_n)-\omega^k_nA_2(\omega^{2k}_n)=A_1(\omega^k_{n\over2})-\ \omega^k_nA_2(\omega^k_{n\over2}) =A1?(ωn2k?)?ωnk?A2?(ωn2k?)=A1?(ω2n?k?)?ωnk?A2?(ω2n? k?)

• What was found?

A(ωkn)A(ωnk) A(\omega^k_n) A(ωnk?) and A(ωk + n2n)A(ωnk + n2) A(\omega^{k + {n\over2}} _n) A(ωnk + 2n) Two polynomialsThe only difference in the last thing is the sign
That is, if it is known that A1(ωkn2)A1(ωn2k) A_1(\omega^k_{n\over 2}) A1?(ω2n?k?) and A2(ωkn2)A2(ωn2k) A_2(\ omega^k_{n\over 2}) A2?(ω2n?k?), we can at the same time know A(ωkn)A(ωnk) A(\omega^k_n) The value of A(ωnk?) and A(ωk + n2n)A(ωnk + n2) A(\omega^{k + {n\over2}}_n) A(ωnk + 2n)
Now we can recursively divide and conquer to do FFTFFT FFT FFT

Each time you backtrack, only scan the previous half of the sequence to get the answer to the next half of the sequence.
n==1n==1 n==1 When n==1 there is only one constant term, just returnreturn return return
Time complexityO(nlog2n)O(nlog?2n) O(n\log_2n) O(nlog2?n)

IFFT (Inverse Fast Fourier Transform)

Think about it, we not only need to know FFTFFT FFT FFT, but also IFFT (Inverse Fast Fourier Transform)
We multiply two polynomials (also called convolution). After doing FFTFFT FFT FFT twice, we also know the point value representation of the product polynomial.
But we usually use coefficients to represent polynomials, and point value representation is meaningless.

• How to convert a polynomial in point value representationquickly back to coefficient representation?

• IDFTIDFT IDFT IDFT violent O(n2)O(n2) O(n^2) O(n2) do? In fact, you can also use FFTFFT FFT. FFT takes O(nlog2n)O(nlog?2n) O(n\log_2n) O(nlog2?n) time.

Have you ever wondered why Fourier substituted ωknωnk \omega^k_n ωnk? as xx x x instead of some other number?
There are reasons but I can’t understand them even if they are
Since I am a fool, I only have to remember one conclusion

• A polynomial is multiplied by the complex conjugate of the unit root during the divide and conquer process. After dividing and conquering, each term divided by nn n n is the coefficient of each term of the original polynomial

This means that FFTFFT FFT FFT and IFFTIFFT IFFT IFFT canbe done together

Simple version of FFT board

c++c++c++c++ has its own complex template complexcomplex complex complex library
a.real()a.real() a.real() a.real() represents the real part of the complex number aa a a

#include<complex>
#define cp complex<double>

void fft(cp *a,int n,int inv)//inv is the symbol for taking conjugate complex numbers
{
    if (n==1)return;
    int mid=n/2;
    static cp b[MAXN];
    fo(i,0,mid-1)b[i]=a[i*2],b[i + mid]=a[i*2 + 1];
    fo(i,0,n-1)a[i]=b[i];
    fft(a,mid,inv),fft(a + mid,mid,inv);//Divide and conquer
    fo(i,0,mid-1)
    {
        cp x(cos(2*pi*i/n),inv*sin(2*pi*i/n));//inv depends on whether it is a conjugate complex number
        b[i]=a[i] + x*a[i + mid],b[i + mid]=a[i]-x*a[i + mid];
    }
    fo(i,0,n-1)a[i]=b[i];
}

Multiply two polynomials a,ba,b a,b a,b and then convert the coefficient polynomial cc c c. Usually you only need to type this short paragraph

 cp a[MAXN],b[MAXN];
int c[MAXN];
fft(a,n,1),fft(b,n,1);//1 coefficient turning point value
fo(i,0,n-1)a[i]*=b[i];
fft(a,n,-1);//-1 point value to coefficient
fo(i,0,n-1)c[i]=(int)(a[i].real()/n + 0.5);//Pay attention to accuracy

It is obvious thatFFT can only handle polynomialsnn n nthat are integer powers of22 2 2
So before FFTFFT FFT FFT be sure to adjust nn n n to the power of 22 2 2

fo(i,0,n-1)rev[i]=(rev[i>>1]>>1)|((i & amp;1)<<(bit-1));

void fft(cp *a,int n,int inv)
{
    int bit=0;
    while ((1<<bit)<n)bit + + ;
    fo(i,0,n-1)
    {
        rev[i]=(rev[i>>1]>>1)|((i & amp;1)<<(bit-1));
        if (i<rev[i])swap(a[i],a[rev[i]]);//If you don’t add this if, it will be exchanged twice (that is, there is no exchange)
    }
    for (int mid=1;mid<n;mid*=2)//mid is half the length of the sequence to be merged
    {
    cp temp(cos(pi/mid),inv*sin(pi/mid));//The unit root, the coefficient 2 of pi has been eliminated
        for (int i=0;i<n;i + =mid*2)//mid*2 is the length of the sequence to be merged, i is the bit to be merged into
{
            cp omega(1,0);
            for (int j=0;j<mid;j + + ,omega*=temp)//Scan only the left half and get the answer to the right half
            {
                cp x=a[i + j],y=omega*a[i + j + mid];
                a[i + j]=x + y,a[i + j + mid]=x-y;//This is the butterfly transformation or something
            }
        }
    }
}