Fast Fourier transform (FFT) (next) __ Fast Fourier transform

About learning the data of FFT algorithm the most recommended or the 30th chapter of Algorithmic Introduction (third edition), the polynomial and fast Fourier transform, the basic knowledge is very comprehensive.

Basic concept of FFT algorithm:

FFT (Fast Fourier transformation), which is rapid Fourier transform, is an accelerated algorithm for discrete Fourier transform, which can complete the DFT in O (Nlogn) time, and can complete the inverse DFT in the same complexity time by using the similarity.
DFT (discrete Fourier Transform) is a discrete Fourier transform, where the coefficients vector of the polynomial is converted to a point-value representation process.
In the ACM-ICPC competition, FFT algorithms are often used to accelerate polynomial multiplication, that is, in O (nlogn) complexity to complete the polynomial multiplication, of course, the actual application is not limited to these, often need to construct polynomial multiplication to count the problem, Also need to use FFT algorithm to solve , some related problems in this article will also mention.

Basic mathematical knowledge required by FFT algorithm:
Polynomial-Related:
Polynomial-related definitions: a polynomial with x as a variable is defined on an algebraic field F, and function A (x) is expressed as a form and:

Called coefficients, all coefficients belong to algebraic domain F, if the maximum non-0 coefficient of a polynomial is, then the number of times this polynomial is K, degree (A) =k, any integer that is strictly greater than the number of polynomials is the number bounds of the polynomial.
On polynomial addition and multiplication believe that the reader of this blog will be the most basic middle school algorithm, the time complexity of polynomial addition is O (n), and the time complexity of multiplication is O (n^2) (n is the number of two polynomials A and b), if the traditional method of calculation is used.

The expression of a polynomial:
In ordinary study, the most common is the expression of the coefficients of polynomials, the number of times bounded by the number of n polynomial is expressed as a vector of coefficients. But the polynomial also has a more commonly used representation method, namely the polynomial point value expression way.

The point-value expression of a polynomial with N in a number bounded is a set of n-point pairs

Makes the k=0,1 of arbitrary integers,.. n−1, each is different, and =a ().

The correctness proof of the point value expression:
For any set of N-point-value pairs, if there is a polynomial a (x) with a number bounded by N, then the n points, then:

The n*n matrix of the leftmost one is called the Vandermonde matrix, and it can be proved by mathematical induction that its determinant value is:.

When the value of the determinant is not 0 when 22 is different, the matrix is reversible, so there exists a unique solution, so the point value expression of the polynomial is reasonable.
The corresponding method of determining the values of polynomial coefficients directly through the coordinates of N-Points is present, interested readers can query the relevant data of Lagrange formula, using Lagrange interpolation formula can get polynomial expression in the time complexity of O (n^2), this is also a problem in the introduction of algorithm

The multiplication of polynomials under point-value expressions is not difficult to find, if the point values of polynomial a (x), B (x) are expressed respectively:

And

So if the polynomial C (x) =a (x) B (x), then the point value of C (x) is expressed as:

convolution :
For the coefficients vector of the two polynomials and the coefficients vector of the polynomial multiplied by the two polynomial, the coefficient vector c is the convolution of the input vectors A and B, and is recorded as C=a⊗b.

In the simple method of polynomial multiplication, the coefficients of each polynomial are calculated by means of convolution under the representation of coefficients. The time complexity is O (n^2), but the FFT is to first convert the polynomial from the representation of the coefficients to the point-value notation (which can be done in O (Nlogn) time complexity, That is, the accelerated DFT transformation, then the product calculation under the point-value notation, the point-value representation of the result is obtained in the time complexity of O (n), then the inverse DFT transformation is performed, and the inverse DfT transform is obtained by the N*logn of O (the time complexity) .

To understand the DFT, you need a certain number of mathematical knowledge:

Complex numbers related to the basics:

unit plural root : n the plural root of a unit is the satisfying of all complex numbers ω, n times the number of units complex root just have n, they are, where I is the plural unit, k=0,1,2...n−1, in the complex plane of the n root evenly distributed on the radius of 1 circle, the definition of the complex number index is as follows:

Which is called the primary N-Time unit root (this definition seems to be useless next)
On several theorems and lemma of complex roots:

Elimination lemma : for any integer there

Prove:

A corollary: for arbitrary even n > 0 have

Proof: Set n = 2*k so

binary lemma : If n > 0 is an even number, then the set of the square of N-N-units complex roots is the set of n/2 of the plural roots of N/2 units

Proof: In fact, this lemma is proof

The binary lemma has a great effect on the conversion of polynomial coefficients to point values by the use of partition, which guarantees that the recursive sub problem is half of the original problem scale.

summation lemma : nonnegative integer k n≥1 to arbitrary integers and not divisible by N, having

This problem can be obtained by geometric progression sum formula:

DFT and FFT, and inverse DFT:

Dft:
In the DfT transformation, we want to calculate the value of polynomial a (x) at the complex root, that is to ask

The vector y= (y0,y1,..., yn−1) is a discrete Fourier transform of the coefficients vector a= (a0,a1,..., an−1), which is recorded as

Fft:

The complexity of the direct computing DfT is O (n^2), the use of the special properties of the complex root can be completed in O (N*logn) time, this method is the FFT method, in the FFT method using a partition strategy to operate, mainly using the elimination of the lemma after the inference.

In the FFT strategy, the number of polynomials is 2 of the integer power, the shortfall is preceded by 0, each step will be the current polynomial a (x), the number of times is 2 multiples, divided into two parts:


And then there's the

So if we can find out the polynomial of the number bounds and the value of the square of the N/2 of n n times, that is, in

, then according to the binary lemma, the number of n is actually only n/2 a different value, that is to say, for each of the two number of times bounded n/2 polynomial, just ask for its N/2 a different value, then the problem is recursive to the original size of half, that is, if you know the results of two child problems, The current problem can be resolved within the complexity of the sum of two sub problems, so the complexity of the recursive problem will be O (Nlogn), with a= (A0,A1,..., an−1) representing the vector of coefficients, y= (y0,y1,..., yn−1) to represent the vector after the discrete transformation, Here gives the C + + code to translate the Code on the Guide (to solve the algorithm in the third edition of the 30th chapter of a problem, to find (0, 1, 2, 3) of the DFT as an example):

#include <bits/stdc++.h> using namespace std;
Const double EPS (1E-8);
 
typedef long Long Lint;
Const double PI = ACOs (-1.0);
    * * This is a recursive implementation of the FFT test, test exercises in the DFT (0, 1, 2, 3)/struct Complex {double real, image;
        Complex (double _real, double _image) {real = _real;
    image = _image;
 
} Complex () {}}; Complex operator + (const Complex &AMP;C1, const Complex &AMP;C2) {return Complex (c1.real + c2.real, C1.image + c2.im
Age); } Complex Operator-(const Complex &AMP;C1, const Complex &AMP;C2) {return Complex (c1.real-c2.real, c1.image-c
2.image); } Complex operator * (const Complex &AMP;C1, const Complex &AMP;C2) {return Complex (c1.real*c2.real-c1.image*c2.i
Mage, C1.real*c2.image + c1.image*c2.real);
    } complex* Recursivefft (Complex a[], int n)//n represents the dimension of vector a {if (n = 1) return A;
    Complex WN = Complex (cos (2*pi/n), sin (2*pi/n));
    Complex w = Complex (1, 0);
    complex* a0 = new complex[n >> 1]; complex* A1 = new complex[n >> 1];
        for (int i = 0; i < n; i++) if (I & 1) a1[(i-1) >> 1] = a[i];
    else a0[i >> 1] = a[i];
    Complex *y0, *y1;
    y0 = Recursivefft (a0, n >> 1);
    Y1 = recursivefft (a1, n >> 1);
    Complex* y = new Complex[n];
        for (int k = 0; k < (n >> 1); k++) {y[k] = Y0[k] + w*y1[k];
        Y[k + (n >> 1)] = Y0[k]-w*y1[k];
    W = w*wn;
} return y;
    int main () {complex* a = new complex[10];
    A[0] = Complex (0, 0);
    A[1] = Complex (1, 0);
    A[2] = Complex (2, 0);
    A[3] = Complex (3, 0);
    complex* ans = new complex[10];
    Ans = Recursivefft (A, 4); for (int i = 0; i < 4 i++) cout<< "(" <<ans[i].real<< ")" << "+" << "(" <<ans[i].im
    
    age<< ")" << "I" <<endl;
return 0; }

The coefficients vectors (0, 1, 2, 3) can be obtained after the DFT (6, -2-2i,-2, -2+2i)

What is needed to be noted in this line of sight is the y[k] = Y0[k] + w*y1[k] and y[k + (n >> 1)] = Y0[k]-w*y1[k]; the change is taking advantage of the even number of N, having been established, to get the array y of the Recursivefft All values

There is also a conceptual thing: and all set up, in which it will be called the rotational factor

After the DfT operation, the expression of the coefficients of polynomial is transformed to the point value, then how to complete the inverse DFT in the time of O (Nlogn), the point value expression is reduced to the coefficient expression by inverse DFT.

Inverse DFT:

According to the relationship between vector y and coefficient vector A of DFT, we can express the relation between them by the form of matrix product, i.e.

To restore the vector y of the DfT change to vector A, we need to use the inverse matrix to multiply the vector y, which requires a theorem.

In the matrix, it is not difficult to find for arbitrary, then can find such a matrix, for arbitrary, (J, K) out of the elements, is the matrix of the inverse matrix.

The proof is as follows: to prove the reciprocal of these two matrices, it is proved that the product is a unit matrix, consider the product of two matrices in (J, j′) of the elements, you can find that this element is

When j′=j, this sum is 1, otherwise, according to the summation lemma, this sum is 0, so the two matrices are mutually inverse

Then, based on this inverse matrix, we can find out what to compute, there is a relational comparison between this formula and the previous DFT inside Y and a, you can find that only need to replace it, the final result needs to be divided by N, so the method of calculating the inverse DFT and computing DFT and similar, can be in O (NLOGN) Time complexity within the solution.

Convolution theorem:
For any two vectors a and b of length n, where n is the power of 2, there

Where vectors A and b are filled with 0, so that the length reaches 2n, and the dot multiplication (that is, the number on each dimension) that ⋅ represents two 2n elements is multiplied.
This equation is actually the result of the convolution operation of the polynomial in the multiplication, which is equivalent to the process of multiplying the coefficients by the DfT transformation into the point-value expression and then swapping back.

Iterative implementation of FFT algorithm:

In the FFT algorithm for recursive DFT, we use the binary lemma to reduce the size of the computation each time by dividing the coefficients vector a into two parts, and we can find that in each grouping they satisfy the grouping of such a complete binary tree (n is the Power of 2):

Recursive coefficients grouping of FFT

By the flow of the diagram above, it can be seen that the order of the sub nodes of the last layer is actually the order in which the subscript is converted into a binary string in a dictionary sequence.

For example, A0~a7 gets the sequence of A0, A4, A2, A6, A1, A5, A3, A7 subscript binary is 000, 100, 010, 110, 001, 101, 011, 111 the corresponding string of the reverse is 000, 001, 010, 011, 100, 10 1, 110, 111 This reverse binary is just 0, 1, 2, 3, the binary representation of 4, 5, 6, 7, so that we can get the subscript order of the following layer within the complexity of O (Nlogn), and then iterate up the values of the father's node based on the result of the child nodes, so that the calculation will directly avoid Recursion, if the direct recursion in some OJ may cause a stack of errors, so it is better to use an iterative approach.

On the iterative form of the FFT algorithm, C + + code implementation is as follows (also for (0, 1, 2, 3) of the DFT transformation as an example):

This code also carries out a reverse DFT, similar to the process of DFT and inverse DFT, plus a marker to determine which one is currently executing.

#include <bits/stdc++.h> using namespace std;
Const double EPS (1E-8);
 
typedef long Long Lint;
 
Const double PI = ACOs (-1.0);
    struct Complex {double real, image;
        Complex (double _real, double _image) {real = _real;
    image = _image;
 
} Complex () {}}; Complex operator + (const Complex &AMP;C1, const Complex &AMP;C2) {return Complex (c1.real + c2.real, C1.image + c2.im
Age); } Complex Operator-(const Complex &AMP;C1, const Complex &AMP;C2) {return Complex (c1.real-c2.real, c1.image-c
2.image); } Complex operator * (const Complex &AMP;C1, const Complex &AMP;C2) {return Complex (c1.real*c2.real-c1.image*c2.i
Mage, C1.real*c2.image + c1.image*c2.real);
    int rev (int id, int len) {int ret = 0; for (int i = 0; (1 << i) < Len;
        i++) {ret <<= 1;
    if (ID & (1 << i)) RET |= 1;
return ret; //When dft= 1 o'clock is a DFT, the DFT =-1 is the inverse DFT complex* iterativefft (complex* A, int len, int DFT//DFT transform for an array of Len (2 power) {complex* a = new complex[len];//after storing an array of array A in the order of a group (int i = 0; i < len; i++)
    A[rev (i, len)] = A[i]; for (int s = 1; (1 << s) <= Len;
        s++) {int m = (1 << s);
        Complex wm = Complex (cos (dft*2*pi/m), sin (dft*2*pi/m));
            for (int k = 0; k < Len; k + m) This layer contains the number of array elements (1 << s) {Complex w = Complex (1, 0);  for (int j = 0; J < (M >> 1); j +)/binary lemma, compute father Node {Complex T = w*a[k + J based on two child nodes
                + (M >> 1)];
                Complex u = a[k + j];
                A[k + J] = U + t;
                A[k + j + (M >> 1)] = u-t;
            W = W*WM;
    }} if (DFT = 1) for (int i = 0; i < len; i++) a[i].real/= len, a[i].image/= Len;
return A;
    int main () {complex* a = new complex[4];
    A[0] = Complex (0, 0);
    A[1] = Complex (1, 0);
    A[2] = Complex (2, 0); A[3] = Complex (3, 0);
    A = Iterativefft (A, 4, 1);
    cout<< "----------after DFT----------" <<endl;
    for (int i = 0; i < 4; i++) printf ("%.9f + (%.9f) i\n", A[i].real, A[i].image);
    cout<< "----------after DFT-1----------" <<endl;
    A = Iterativefft (A, 4,-1);
    for (int i = 0; i < 4; i++) printf ("%.9f + (%.9f) i\n", A[i].real, A[i].image);
return 0; }

Through the example of the above code, the implementation of the FFT algorithm is basically no problem, in addition to the introduction of algorithm some of the problems are very good, easy to familiarize yourself with the many details of this algorithm, here is not mentioned.
After these studies, it is necessary to carry out some practical exercises, followed by the relevant exercises section.

HDU 1402 a*b Problem Plus Large integer Multiplication
HDU 4609 3-idiots FFT count
Uvalive 4671 K-neighbor substrings FFT calculation string Hamming distance
UVA 12298 Super Poker II FFT count, long double
URAL 1996 Cipher Massage 3 FFT + KMP
Codechef Countari Arithmetic progressions FFT + chunking
Zoj 3856 Goldbach FFT count
Uvalive 6886 Golf Bot FFT template problem
Hdu 4093 Xavier is Learning to Count tolerance principle + FFT (2011 Shanghai Live game c)

HDU 5751 bestcoder Round #84 eades (segment tree +fft)

HDU 5730 2016 multi-school 1 Shell Necklace (CDQ Division +fft)

2016 ACM Hong Kong network Race A question A+b Problem (FFT)

Later there are topics can see my blog, may have. Finally, thank Ichimei for the explanation of FFT. Orz ..... It really helps to understand the FFT.