FFT learning notes (I think it is detailed)

introduce

What is? \(\text{FFT}\) ？
Anyway, I see \(\text{wiki}\) There are a lot of strange things on it .

The fast Fourier transform （ English ：Fast Fourier Transform, FFT）, It is the discrete Fourier transform of fast computing sequence （DFT） Or its inverse transformation . Fourier analysis transforms the signal from the original domain （ Usually time or space ） Convert to the representation of the frequency domain or reverse .FFT Will pass DFT The matrix is decomposed into sparse （ Mostly zero ） Product of factors to quickly calculate such transformations .—— \(\text{wikipedia}\)

Anyway, I have no brain, I can't understand .
For me, ,\(\text{FFT}\) Is to be able to multiply polynomials from \(O(n^2)\) become \(O(n\log n)\) Fairy stuff .

Text

Coefficient representation and point value representation

For the coefficient representation , Is to use the coefficients of a polynomial to express this polynomial .
for instance ：

\[f(x)=a_1x^3+a_2x^2+a_3x+a_4\Leftrightarrow f(x)=\{a_1,a_2,a_3,a_4\} \]

So for the point value representation , Correspondingly, the polynomial is represented by several points on the function .
Students who have studied primary school mathematics must know ：\(n+1\) One point determines one \(n\) Sub polynomial .
If you prove it, you can consider mathematical induction ./xyx
Take another example , The point value representation is like this ：

\[f(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n⇔f(x)=\{(x_0,y_0),(x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)\}\\ \]

As mentioned above, it is necessary to convert the coefficient representation into the point value representation . So why is this ？
Now let's show the polynomial multiplication of point-valued representation ：

\[f(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),\cdots,(x_n,f(x_n))\}\\ g(x)=\{(x_0,g(x_0)),(x_1,g(x_1)),(x_2,g(x_2)),\cdots,(x_n,g(x_n))\}\\ F(x) = f(x)\times g(x)\\ F(x)=\{(x_0,g(x_0)\times f(x_0)),(x_1,g(x_1)\times f(x_1)),\cdots,(x_n,g(x_n)\times f(x_n))\} \]

The plural

\( The plural = The set of real Numbers + imaginary number \)
Be honest , Direct delivery of dry goods , We define :

\[i=\sqrt{-1} \]

In this way, we can express numbers that we cannot express in the range of real numbers .
So how to express a plural number ：

\[Num=a+bi\ \ (a,b \in \mathbb{R}) \]

Then we put \(Num=a+bi\) As a function , hold \(a\) and \(b\) They correspond to each other \(x\) Axis and \(y\) Axis .
You can get the complex plane , It's like this ：

The abscissa is the real number axis , The ordinate is the imaginary axis , In this way, each imaginary number can be regarded as a vector .
Corresponding , Imaginary numbers can be expressed in ordinary coordinates and polar coordinates ：

\[(x,y)\quad and \quad (r,\theta) \]

The meaning of multiplying two complex numbers is given below ：

\[\begin{split}&\quad (a+bi)\times(c+di)\\&= ac+adi+bci+bdi^2\\ &=(ac+bd)+(ad+bc)i \end{split} \\\quad \\ (r_1,\theta_1)\times(r_2,\theta_2)=(r_1\times r_2,\theta_1+\theta_2) \]

\(\tt DFT\) （ Discrete Fourier transform ）

Now it's over Point value representation and The plural Knowledge about , Next is the dry goods part .

Above, we have illustrated the convenience of point-valued representation in polynomial multiplication through such an example .
Next, let's look at how to convert polynomials from coefficient representation to point value representation , This process is called \(\text{DFT}\) .

The so-called point value representation , That is to say \(n\) Take on the polynomial \(n+1\) A little bit , To represent .
Formalized , It can be expressed as ：

\[F(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}+a_nx^n\\ \rightarrow F(x)=\{(x_0,F(x_0)),(x_1,F(x_1)),(x_2,F(x_2)),\cdots,(x_n,F(x_n))\} \]

Then you can be surprised to find , Take a few with you \(x_i\) Go in and calculate \(F(x_i)\) Just fine .
But if you graduate from primary school , You can find that it's better to be direct \(O(n^2)\) violence .

So what should I do ？
We wonder if there are some \(x\) bring \(x^n\ (n\in \tt Z^+)\) The result is \(1\) .
This seems to be a very good idea , But how many such numbers are there ？
I can blurt out two \(1\) and \(-1\) , Think about it and you will find that \(i\) and \(-i\) Can also be .

But after careful consideration （ Look at the explanation ） It can be found that all points on the unit circle in the figure below meet the conditions .

For convenience , We are taking this \(n\) This unit will be divided equally at a point .

We from \((1, 0)\) This point starts , From... In a counterclockwise direction \(0\) Start numbering , Form like \(\omega_n^k\) .
among \(n\) It means that a total of \(n\) A little bit ,\(k\) Represents the number of the current point .

By the complex multiplication we introduced before Multiplication of module length , Add degrees ：

\[(r_1,\theta_1)\times(r_2,\theta_2)=(r_1\times r_2,\theta_1+\theta_2) \]

And combine the properties of the unit circle （ The distance from all points to the origin is \(1\)）.
You can get from \(\omega_n^1\) The switch to \(\omega_n^k\) Formula ：

\[(\omega_n^1)^k=\omega_n^k \]

We call \(\omega_n^1\) by \(n\) Subunit root .
So we can find that , Let's go straight into \(\omega_n^i\) That's all right. .

Some useful properties of unit root

Before we know the nature of everything , We need to know the unit root first \(\omega_n^i\) How to say ：

\[\omega_n^k=\cos\frac{k}{n}2\pi+i\times\sin\frac{k}{n}2\pi \]

For the proof of this thing, you can draw a point directly on the unit circle, and then you can get a basic knowledge of trigonometric functions .

Property one

\[\omega_n^k=\omega_{2n}^{2k} \]

If you prove it, just follow the formula set given above , And then I found that I can cut scores .
Then I think we can get further ：

\[\omega_n^k=\omega_{Pn}^{Pk} \quad (p\in \tt Z^+) \]

Obviously, but it doesn't seem to be of much use .

Property two

\[\omega_n^{k+\frac{n}{2}}=-\omega_n^k \]

Write down the proof a little ：

\[\omega_n^k=\cos\frac{k}{n}2\pi+i\times\sin\frac{k}{n}2\pi\\ \begin{split} \omega_n^{k+\frac{n}{2}}&= \cos\frac{k+\frac{n}{2}}{n}2\pi+i\times\sin\frac{k+\frac{n}{2}}{n}2\pi\\ &=\cos(\frac{k}{n}2\pi+\pi)+i\times\sin(\frac{k}{n}2\pi+\pi) \end{split} \]

If it turns into this step, we won't prove it in the next step , I still don't understand the suggestion to retake junior high school mathematics .

Property three

\[\omega_n^0=\omega_n^n \]

Relatively simple , I won't say why .

\(\tt FFT\) （ The fast Fourier transform ）

Here he comes , Here he comes , Until now, he finally came ....

As mentioned before, we can directly bring \(\omega_n^i\) To calculate the point value .
Yes , I think this method is efficient , clever , Doggerel , It embodies human wisdom .
But wait. , Although the process of calculating coefficients is eliminated , But for each \(\omega_n^i\) We still have to \(O(n)\) Calculate the result .
Then I moved around and counted my fingers , I found that there were \(n\) individual \(\omega_n^i\) Value , And then \(O(n^2)\) 了 .
So what should we do ？

Look at the solution carefully , It can be found from the perspective of divide and conquer .
Be careful ： The following guarantees \(n\) by \(2\) Integer power of .
Let's set a polynomial ：

\[\begin{split} F(x)&=\sum_{i=0}^{n-1}a_ix^i\\ &=a_0 + a_1x + a_2x^2 + \cdots + a_{n-1}x^{n-1} \end{split} \]

Then find a way to \(F(x)\) Divided into two parts .
The method adopted here is based on \(F(x)\) The parity of the subscript is divided into two parts .

\[\begin{split} F(x) &=a_0 + a_1x + a_2x^2 + \cdots + a_{n-1}x^{n-1}\\ &=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots+a_{n-1}x^{n-2}) \end{split} \]

Next, we find that the structure of these two polynomials is exactly the same .
Let's set these two polynomials as \(F_1(x)\) and \(F_2(x)\) .

\[F_1(x)=a_0+a_2x^2+\cdots+a_{n-2}x^{n-2}\\ F_2(x)=a_1+a_3x^2+\cdots+a_{n-1}x^{n-2}\\ F(x)=F_1(x)+xF_2(x) \]

It is found that such coefficients are discontinuous , Not so perfect , So we changed again .

\[F_1(x)=a_0+a_2x^1+\cdots+a_{n-2}x^{\frac{n}{2}-1}\\ F_2(x)=a_1+a_3x^1+\cdots+a_{n-1}x^{\frac{n}{2}-1}\\ F(x) = F_1(x^2)+xF_2(x^2) \]

At this time, you can find that this form is very beautiful .
The next step is to bring it directly \(\omega_n^i\) Operating the .
We then set \(k<\frac{n}{2}\) And then put \(\omega_n^k\) Directly into .

\[\begin{split} F(\omega_n^k) &= F_1((\omega_n^k)^2)+\omega_n^kF_2((\omega_n^k)^2)\\ &=F_1(\omega_{n}^{2k})+\omega_n^kF_2(\omega_{n}^{2k})\\ &=F_1(\omega_{\frac{n}{2}}^{k})+\omega_n^kF_2(\omega_{\frac{n}{2}}^{k}) \end{split} \]

The first step is to directly bring , If there is a problem, the primary school suggests to rebuild .
I wrote the second step before , The formula is like this ：

\[(\omega_n^1)^k=\omega_n^k \]

Of course , The application here is universal , If you have any questions, just push .
As for the third step , I think it will be faster to calculate the proportion directly .
about \(F(\omega_n^{k+\frac{n}{2}})\) Directly into ：

\[\begin{split} F(\omega_n^{k+\frac{n}{2}}) &= F_1((\omega_n^{k+\frac{n}{2}})^2)+\omega_n^{k+\frac{n}{2}}F_2((\omega_n^{k+\frac{n}{2}})^2)\\ &=F_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}F_2(\omega_n^{2k+n})\\ &=F_1(\omega_n^{2k}\omega_n^n)-\omega_n^{k}F_2(\omega_n^{2k}\omega_n^n)\\ &=F_1(\omega_n^{2k})-\omega_n^{k}F_2(\omega_n^{2k})\\ &=F_1(\omega_{\frac{n}{2}}^{k})-\omega_{n}^{k}F_2(\omega_{\frac{n}{2}}^{k})\\ \end{split} \]

It is troublesome to introduce each step one by one , Let's just hand it over or turn to the previous formula .
Observe the first formula and the second formula , The only difference I found was the symbol .
Then divide and conquer directly to solve , Time complexity \(O(n\log n)\) .

\(\tt IFF\) ( Fast Fourier inverse transform )

That is to convert the point value representation into the coefficient representation we want .
Here's the conclusion , It's disgusting to prove that it's true , So I won't prove it .

A polynomial is multiplied by the conjugate complex of the unit root in the process of divide and conquer , Divide and conquer each item by \(n\) That is, the coefficients of each term of the original polynomial

That is to do it again \(\tt FFT\) Divide each bit by \(n\) That's all right. .

Code implementation and optimization

Code Complex type encapsulation

struct cp {
  double x, y;
  cp (double xx = 0, double yy = 0) {x = xx; y = yy;};

  friend cp operator +(cp p, cp q) {return cp(p.x + q.x, p.y + q.y);}
  friend cp operator -(cp p, cp q) {return cp(p.x - q.x, p.y - q.y);}
  friend cp operator *(cp p, cp q) {return cp(p.x * q.x - p.y * q.y, p.y * q.x + p.x * q.y);}
}a[N], b[N];

Code No optimization

I didn't write the code , Anyway, it is simulated according to the previous formula , Just have a look .

Click to view the code

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

void fft(cp *a, int n, int inv) //inv Is the sign of the conjugate complex 
{
  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 the conjugate complex number is taken 
    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];
}

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 conversion coefficient 
fo(i, 0, n - 1)c[i] = (int)(a[i].real() / n + 0.5); // Pay attention to precision 

Be careful ：\(\tt FFT\) Before you do that, you have to \(n\) Modulation \(2\) Omega to an integer power .
Obviously, the one above can't even pass the template question .
So here we need to consider how to optimize \(\tt FFT\) .

Look at the original sequence and the inverted sequence , The sequence to be solved is actually the binary inversion of the subscript of the original sequence ！
Therefore, the process of classifying sequences according to the parity of subscripts is actually unnecessary .
So we can \(O(n)\) Using some kind of operation to get the sequence we require , Then continue to merge upward .
—— \(\tt luogu\) Solution to a problem

Code There is optimization , Passable

Click to view the code

#include <bits/stdc++.h>

#define file(a) freopen(a".in", "r", stdin), freopen(a".out", "w", stdout)

#define Enter putchar('\n')
#define quad putchar(' ')

#define N 3000005

namespace IO {

template <class T>
inline void read(T &a);
template <class T, class ...rest>
inline void read(T &a, rest &...x);

template <class T>
inline void write(T x);

}

struct cp {
  double x, y;
  cp (double xx = 0, double yy = 0) {x = xx; y = yy;};

  friend cp operator +(cp p, cp q) {return cp(p.x + q.x, p.y + q.y);}
  friend cp operator -(cp p, cp q) {return cp(p.x - q.x, p.y - q.y);}
  friend cp operator *(cp p, cp q) {return cp(p.x * q.x - p.y * q.y, p.y * q.x + p.x * q.y);}
}a[N], b[N];

const double pi = acos(-1.0);

int n1, n2, n, rev[N], c[N];

inline void FFT(cp *a, int n, int inv) {
  int bit = 0;
  while ((1 << bit) < n) bit++;
  for (int i = 1; i < n; ++i) {
    rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
    if (i < rev[i])
      std::swap(a[rev[i]], a[i]);
  }
  for (int mid = 1; mid < n; mid <<= 1) {
    cp temp(cos(pi / mid), inv * sin(pi / mid));
    for (int i = 0; i < n; i += mid * 2) {
      cp omega(1,0);
      for (int j = 0; j < mid; ++j, omega = omega * temp) {
        cp x = a[i + j], y = omega * a[i + j + mid];
        a[i + j] = x + y;
        a[i + j + mid] = x - y;
      }
    }
  }
}

signed main(void) {
  // file("P3803");
  IO::read(n1, n2);
  n = std::max(n1, n2);
  for (int i = 0, num; i <= n1; ++i) IO::read(num), a[i].x = num;
  for (int i = 0, num; i <= n2; ++i) IO::read(num), b[i].x = num;
  n = n1 + n2;
  for (int i = 0; i <= 30; ++i)
    if ((1 << i) > n) {
      n = (1 << i);
      break;
    }

  FFT(a, n, 1); FFT(b, n, 1);
  for (int i = 0; i < n; ++i) a[i] = a[i] * b[i];
  FFT(a, n, -1);
  for (int i = 0; i <= n1 + n2; ++i)
    c[i] = (int)(a[i].x / n + 0.5);
  
  for (int i = 0; i <= n1 + n2; ++i)
    IO::write(c[i]), quad;
  Enter;
}

namespace IO {

template <class T>
inline void read(T &a) {
  T s = 0, t = 1;
  char c = getchar();
  while ((c < '0' || c > '9') && c != '-')
    c = getchar();
  if (c == '-')
    c = getchar(), t = -1;
  while (c >= '0' && c <= '9')
    s = (s << 1) + (s << 3) + (c ^ 48), c = getchar();
  a = s * t;
}
template <class T, class ...rest>
inline void read(T &a, rest &...x) {
  read(a);
  read(x...);
}

template <class T>
inline void write(T x) {
  if (x == 0) putchar('0');
  if (x < 0) putchar('-'), x = -x;
  int top = 0, sta[55] = {0};
  while (x) 
    sta[++top] = x % 10, x /= 10;
  while (top)
    putchar(sta[top] + '0'), top--;
  return ;
}

}

Recommend here Some Zhihu column , hold \(\tt FFT\) Optimization is very clear .

\(\tt NTT\) I'll still see it , however \(\tt FFT\) Give me some collapse ...