Fast Fourier Transform (FFT) Essays

Finally, I learned FFT and recorded a piece of water with my notes

There are a lot of things on the front-end knowledge network, so let's not go into details here, just get to the point

 

01 Introduction to FFT

This only points out the general application of FFT in competition, that is, optimizing polynomial multiplication

Generally, calculating the result of multiplying two polynomials of size $n$is $O(n^2)$, but the magical FFT optimizes it to $O(nlogn)$

The FFT process is generally:

Coefficient representation of polynomial $\longrightarrow$Point value representation of polynomial $\longrightarrow$Coefficient representation of polynomial

Every step on the Internet is called differently. Here the first step is called Fast Fourier Transform and the second is Fast Fourier Inverse Transform.

 

02 Fast Fourier Transform

First point out that each of the next $n$is an integer power of $2$

First we have a polynomial of $n$term expressed by a known coefficient

$A(x)=a_0+a_1x+a_2x^2+\dots+a_{n-1}x^{n-1}$

To determine its point value expression $(y_0,y_1,y_2,\dots,y_{n-1})$, it is simple to substitute $n$different values, which is obviously $O(n^2)$

The following describes the Fast Fourier Transform approach

First classify polynomials by parity

$A(x)=(a_0+a_2x^2+\dots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+\dots+a_{n-1}x^{n-1})$

$A(x)=(a_0+a_2x^2+\dots+a_{n-2}x^{n-2})+x\cdot(a_1+a_3x^2+\dots+a_{n-1}x^{n-2})$

set up

$A_1(x)=a_0+a_2x+\dots+a_{n-2}x^{\tfrac{n-2}{2}}$

$A_2(x)=a_1+a_3x+\dots+a_{n-1}x^{\tfrac{n-2}{2}}$

  Not hard to find

$A(x)=A_1(x^2)+xA_2(x^2)$

Make $k<\frac{n}{2}$

Put $\omega_{n}^k$substituted

$A(\omega_{n}^k)=A_1(\omega_{n}^{2k})+\omega_{n}^{k}A_2(\omega_{n}^{2k})$

$A(\omega_{n}^k)=A_1(\omega_{\tfrac{n}{2}}^{k})+\omega_{n}^{k}A_2(\omega_{\tfrac{n}{2}}^{k})$

Put $\omega_{n}^{k+\frac{n}{2}$substituted

$A(\omega_{n}^{k}+\tfrac{n}{2})=A_1(\omega_{n}^{2k+n})+\omega_{n}^{k+\tfrac{n}{2}}A_2(\omega_{n}^{2k+n})$

$A(\omega_{n}^{k}+\tfrac{n}{2})=A_1(\omega_{n}^{2k}\cdot\omega_{n}^{n})-\omega_{n}^{k}A_2(\omega_{n}^{2k}\cdot\omega_{n}^{n})$

$A(\omega_{n}^{k}+\tfrac{n}{2})=A_1(\omega_{n}^{2k})-\omega_{n}^{k}A_2(\omega_{n}^{2k})$

$A(\omega_{n}^k)=A_1(\omega_{\tfrac{n}{2}}^{k})-\omega_{n}^{k}A_2(\omega_{\tfrac{n}{2}}^{k})$

Obviously, only the constant term is different between the two formulas

When $k$traverses all the values in $[0,\frac{n}{2}-1]$k+\dfrac{n}{2}$all the values in $[\dfrac{n}{2},n-1]$

Therefore, the size of the problem is reduced by half, and we only need to enumerate $k$i n $[0,\dfrac{n}{2}-1]$so that we can calculate all the values of $A(\omega_{n}^i)\quad(i\i n[0, n-1])$

If we know $A_1(x),A_2(x)$at $\omega_ {\tfrac{n}{2}}^0, \omega_ {\tfrac{n}{2}}^1, \dots, \omega_ The value of {\tfrac{n}{2}}^{\tfrac{n}{2}-1}$allows you to calculate $A(x)$in the time of $O(n)$using the above two formulas

Ask for $A_1(x),A_2(x)$is exactly a subproblem of $A(x)$and can be solved recursively

 

03 Fast Fourier Inverse Transform

Above we converted the coefficient representation of a polynomial to the point value representation, and here we will study converting the point value representation of a polynomial to the coefficient representation.

Remember that $(a_0,a_1,\dots,a_{n-1})$is a coefficient vector of $A(x)$and we know that the point value of $A(x)$is expressed as $(A(x_0),A(x_1),\dots,A(x_{n-1})$

Set the vector $(d_0,d_1,\dots,d_{n-1})$to be a coefficient vector, expressed as a point value derived from the Fast Fourier Transform, $(a_0,a_1,\dots,a_{n-1})$

Construct a polynomial $F(x)=d_0+d_1x+d_2x^2+\dots+d_{n-1}x^{n-1}$

Set $(c_0,c_1,\dots,c_{n-1})$to be $F(x)$at $x=\omega_ Point value representation at n^{-k}$i.e. $c_k=F(\omega_n^k)$, that is, $c_k=\sum_{i=0}^{n-1}d_i(\omega_n^{-k})^i$

We know $d_k=A(\omega_n^k)$, that is, $d_k=\sum_{j=0}^{n-1}a_j(\omega_n^k)^j$

Combining the two sum formulas above

$c_k=\sum_{i=0}^{n-1} [\sum_{j=0}^{n-1}a_j(\omega_n^i)^j] (\omega_n^{-k})^i$

$\quad \:=\sum_{i=0}^{n-1} \sum_{j=0}^{n-1}a_j(\omega_n^j)^i (\omega_n^{-k})^i$

$\quad \:=\sum_{j=0}^{n-1} a_j \sum_{i=0}^{n-1} (\omega_n^j \omega_n^{-k})^i$

$\quad \:=\sum_{j=0}^{n-1} a_j \sum_{i=0}^{n-1} (\omega_n^{j-k})^i$

Let's discuss the following sum $\sum_separately {i=0}^{n-1} (\omega_n^{j-k})^i$

$j \neq\ k$

Then the latter sum is converted to an equal ratio sum

$\sum_{i=0}^{n-1} (\omega_n^{j-k})^i=\frac{\omega_n^0 [1-(\omega_n^{j-k})^n]}{1-\omega_n^{j-k}}$

$\qquad \qquad \quad \: \: \:=\frac{1-(\omega_n^{j-k})^n}{1-\omega_n^{j-k}}$

$\qquad \qquad \quad \: \: \:=\frac{1-(\omega_n^n)^{j-k}}{1-\omega_n^{j-k}}$

$\qquad \qquad \quad \: \: \:=\frac{1-1^{j-k}}{1-\omega_n^{j-k}}$

$\qquad \qquad \quad \: \: \:=\frac{0}{1-\omega_n^{j-k}}$

$\qquad \qquad \quad \: \: \:=0$

$j = k$

So $\omega_n^{j-k} = 1$

$\sum_{i=0}^{n-1} (\omega_n^{j-k})^i = n$

From both cases above, we know that the whole formula is only valid if and only if $j = k$and the rest is $0$

So there are

$c_j=a_jn$

$a_j = \frac{c_j}{n}$

So here we have the coefficient expression for $A(x)$

Looking at the whole analysis process, we take the point value of $A(x)$for $(A(x_0),A(x_1),\dots,A(x_{n-1})$as the coefficient representation of a new polynomial, $F(x)$and do a quick Fourier transformation of $F(x)$to get $(c_0,c_1,\dots,c_{n-1})$, then divide by $n to get the coefficient representation of $A(x)$ It should be noted that in fast Fourier transform, $x=\omega_n^k$but $\omega_is substituted in the inverse transformation N^{-k}$

 

04 Implementation

Learn the preceding methods, it will be easy to implement them

For $C(x)=A(x)\cdot B(x)$

Convert both $A(x)$and $B(x)$into point value expressions, $(a_0,a_1,\dots,a_{n-1})$and $(b_0,b_1,\dots,b_{n-1})$

The result is transformed into a coefficient expression of $C(x)$

Post a C++ code for the FFT board Title P3803 on Lough Valley

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#define MAXN 4000006
using namespace std;
class complex
{
public:
    complex(){}
    complex(double a,double b)
    {
        this->a=a;
        this->b=b;
    }
    double a,b;
}a[MAXN],b[MAXN];
complex operator+ (complex x,complex y)
{
    return complex(x.a+y.a,x.b+y.b);
}
complex operator- (complex x,complex y)
{
    return complex(x.a-y.a,x.b-y.b);
}
complex operator* (complex x,complex y)
{
    return complex(x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a);
}
const double pi=acos(-1.0);
void FFT(int l,complex *arr,int f)
{
    if(l==1) return;
    int dl=l>>1;
    complex a1[dl],a2[dl];
    for(int i=0;i<l;i+=2)
    {
        a1[i>>1]=arr[i];
        a2[i>>1]=arr[i+1];
    }
    FFT(dl,a1,f);
    FFT(dl,a2,f);
    complex wn=complex(cos(2.0*pi/l),sin(2.0*pi/l)*f),w=complex(1.0,0.0);
    for(int i=0;i<dl;i++,w=w*wn)
    {
        arr[i]=a1[i]+w*a2[i];
        arr[i+dl]=a1[i]-w*a2[i];
    }
}
int n,m,N;
int main()
{
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++)
        scanf("%lf",&a[i].a);
    for(int i=0;i<=m;i++)
        scanf("%lf",&b[i].a);
    N=1;
    while(N<n+m+1) N<<=1;
    FFT(N,a,1);
    FFT(N,b,1);
    for(int i=0;i<N;i++)
        a[i]=a[i]*b[i];
    FFT(N,a,-1);
    for(int i=0;i<n+m+1;i++)
        printf("%d ",(int)(a[i].a/N+0.5));
    puts("");
    return 0;
}

Nothing to do at leisure. Post another copy of Python's

import numpy as np

pi = np.arccos(-1.0)


def read():
    def get_numbers():
        try:
            read.s = input().split()
            read.s_len = len(read.s)
            if read.s_len == 0:
                get_numbers()
            read.cnt = 0
            return 1
        except:
            return 0

    if not hasattr(read, 'cnt'):
        if not get_numbers():
            return 0
    if read.cnt == read.s_len:
        if not get_numbers():
            return 0
    read.cnt += 1
    return eval(read.s[read.cnt - 1])


n = int(read())
m = int(read())


class Complex:
    # Complex class

    def __init__(self, a=0.0, b=0.0):
        self.a = a
        self.b = b

    def __add__(self, other):
        return Complex(self.a + other.a, self.b + other.b)

    def __sub__(self, other):
        return Complex(self.a - other.a, self.b - other.b)

    def __mul__(self, other):
        return Complex(self.a * other.a - self.b * other.b, self.a * other.b + self.b * other.a)


def fft(num, f, args):
    if num == 1:
        return
    div_num = num >> 1
    a1 = []
    a2 = []
    for i in range(0, num, 2):
        a1.append(args[i])
        a2.append(args[i + 1])
    fft(div_num, f, a1)
    fft(div_num, f, a2)
    wn = Complex(np.cos(2.0 * pi / num), np.sin(2.0 * pi / num) * f)
    w = Complex(1.0, 0.0)

    for i in range(0, div_num):
        args[i] = a1[i] + w * a2[i]
        args[i + div_num] = a1[i] - w * a2[i]
        w = w * wn


aa = []
bb = []
for j in range(0, n + 1):
    aa.append(Complex(float(read()), 0.0))
for j in range(0, m + 1):
    bb.append(Complex(float(read()), 0.0))

nn = 1
while nn < n + m + 1:
    nn <<= 1

for j in range(n + 1, nn):
    aa.append(Complex(0.0, 0.0))
for j in range(m + 1, nn):
    bb.append(Complex(0.0, 0.0))

fft(nn, 1, aa)
fft(nn, 1, bb)

for j in range(0, nn):
    aa[j] = aa[j] * bb[j]
fft(nn, -1, aa)

for j in range(0, n + m + 1):
    print(int(aa[j].a / nn + 0.5), end=' ')

Python is too slow...

 

05 Conclusion

You've learned Fast Fourier Transform at last, and to some extent made up for some of the regrets of the past.

Here's a big man's picture explaining FFT's thinking

 

Here's also a recommendation for Big Guys'blogs

Fast Fourier Transform (FFT) Details - A Pawn in the Wind and Moon - Blog Park (cnblogs.com)

An hour to learn Fast Fourier Transform - Know (zhihu.com)

 

Tags: Math FFT

Posted on Thu, 02 Dec 2021 13:23:00 -0500 by sprocket