Fast Fourier Transform

Prerequisite: Complex Numbers.

This article introduces an algorithm that supports computing the multiplication of two \(n\)-degree polynomials in \(O(n\log n)\) time, which is more efficient than the naive \(O(n^2)\) algorithm. Since the multiplication of two integers can also be treated as polynomial multiplication, this algorithm can also be used to accelerate the computation of large integer multiplication.

Introduction

Let us introduce two polynomials \(A\) and \(B\):

\[ \begin{aligned} A ={}& 5x^2 + 3x + 7 \\ B ={}& 7x^2 + 2x + 1 \\ \end{aligned} \]

The product \(C = A \times B\) of the two polynomials can be solved in \(O(n^2)\) time complexity (where \(n\) is the degree of polynomial \(A\) or \(B\)):

\[ \begin{aligned} C ={}& A \times B \\ ={}& 35x^4 + 31x^3 + 60x^2 + 17x + 7 \end{aligned} \]

Clearly, the coefficient \(c_i\) of polynomial \(C\) satisfies \(c_i = \sum_{j = 0}^i a_j b_{i - j}\). For this naive algorithm, the time complexity for computing each term is \(O(n)\), and there are \(O(n)\) terms, so the overall time complexity is \(O(n^2)\).

Can we accelerate it to reduce its time complexity? If we use the Fast Fourier Transform, we can lower the complexity to \(O(n \log n)\).

Fourier Transform

The Fourier Transform is a method for analyzing signals; it can analyze the components of a signal and also synthesize signals from these components. Many waveforms can serve as signal components, and the Fourier Transform uses sine waves as signal components.

Let \(f(t)\) be a function of time \(t\). Then the Fourier Transform can detect the extent to which a frequency \(\omega\) appears periodically in \(f(t)\):

\[ F(\omega)=\mathbb{F}[f(t)]=\int_{-\infty}^{\infty}f(t)\mathrm{e}^{-\mathrm{i}{\omega}t}dt \]

Its inverse transform is

\[ f(t)=\mathbb{F}^{-1}[F(\omega)]=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)\mathrm{e}^{\mathrm{i}{\omega}t}d\omega \]

The form of the inverse transform is very similar to the forward transform, with the denominator \(2\pi\) exactly matching the period of the exponential function.

The Fourier Transform is equivalent to taking the continuous inner product of the time-domain function with a complex exponential function of period \(2\pi\). The inverse transform is also an inner product.

The Fourier Transform has a corresponding convolution theorem, which can convert convolution in the time domain to multiplication in the frequency domain, and vice versa.

Discrete Fourier Transform

The Discrete Fourier Transform (DFT) is a discrete form of the Fourier Transform in both the time and frequency domains, transforming the time-domain sampling of a signal into its DTFT (discrete-time Fourier transform) frequency-domain sampling.

The Fourier Transform is an integral form of continuous function inner product, while the Discrete Fourier Transform is a summation form of inner product.

Let \(\{x_n\}_{n=0}^{N-1}\) be a sequence satisfying certain finite conditions. Its Discrete Fourier Transform (DFT) is:

\[ X_k=\sum_{n=0}^{N-1}x_n\mathrm{e}^{-\mathrm{i}\frac{2\pi}{N}kn} \]

where \(\mathrm{e}\) is the base of the natural logarithm, and \(i\) is the imaginary unit. This transformation is usually denoted by the symbol \(\mathcal {F}\), i.e.,

\[ \hat{x}=\mathcal{F}x \]

Similar to the integral form, its Inverse Discrete Fourier Transform (IDFT) is:

\[ x_n=\frac{1}{N}\sum_{k=0}^{N-1}X_k\mathrm{e}^{\mathrm{i}\frac{2\pi}{N}kn} \]

It can be written as:

\[ x=\mathcal{F}^{-1}\hat{x} \]

In practice, the normalization coefficients in the DFT and IDFT transformation formulas are not crucial. In the definitions above, the coefficients for DFT and IDFT are \(1\) and \(\frac {1}{N}\), respectively. Sometimes we change both coefficients to \(\frac{1}{{\sqrt{N}}}\).

The Discrete Fourier Transform is still a transformation from the time domain to the frequency domain. Due to the particularity of the summation form, there are other interpretations.

If we treat the sequence \(x_n\) as the coefficients of the polynomial \(f(x)\) for \(x^n\), then the computed \(X_k\) is exactly the point value \(f(\mathrm{e}^{\frac{-2\pi \mathrm{i}k}{N}})\) of the polynomial \(f(x)\) evaluated at the unit root \(\mathrm{e}^{\frac{-2\pi \mathrm{i}k}{N}}\).

This provides another interpretation of the convolution theorem, i.e., performing a special evaluation operation on polynomials. The Discrete Fourier Transform is exactly the evaluation of polynomials at unit roots.

For example, compute:

\[ \dbinom{n}{3}+\dbinom{n}{7}+\dbinom{n}{11}+\dbinom{n}{15}+\ldots \]

Define the function \(f(x)\) as:

\[ f(x)={(1+x)}^n=\dbinom{n}{0}x^0+\dbinom{n}{1}x^1+\dbinom{n}{2}x^2+\dbinom{n}{3}x^3+\ldots \]

Then we can find that evaluating at the fourth unit root \(f(\mathrm{i})\) yields the sequence:

\[ f(\mathrm{i})={(1+\mathrm{i})}^n=\dbinom{n}{0}+\dbinom{n}{1}\mathrm{i}-\dbinom{n}{2}-\dbinom{n}{3}\mathrm{i}+\ldots \]

Thus the following summation exactly cancels out the other terms:

\[ f(1)+\mathrm{i}f(\mathrm{i})-f(-1)-\mathrm{i}f(-\mathrm{i})=4\dbinom{n}{3}+4\dbinom{n}{7}+4\dbinom{n}{11}+4\dbinom{n}{15}+\ldots \]

Therefore, the answer to this math problem is:

\[ \dbinom{n}{3}+\dbinom{n}{7}+\dbinom{n}{11}+\dbinom{n}{15}+\ldots=\frac{2^n+\mathrm{i}(1+\mathrm{i})^n-\mathrm{i}(1-\mathrm{i})^n}{4} \]

This math problem evaluates at unit roots, exactly constituting a Discrete Fourier Transform.

Matrix Formula

Since the Discrete Fourier Transform is a linear operator, it can be described by matrix multiplication. In matrix representation, the Discrete Fourier Transform is expressed as:

\[ \begin{bmatrix} X_{0} \\ X_{1} \\ X_{2} \\ \vdots \\ X_{N-1} \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \alpha & \alpha^{2} & \cdots & \alpha^{N-1} \\ 1 & \alpha^{2} & \alpha^{4} & \cdots & \alpha^{2(N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \alpha^{N-1} & \alpha^{2(N-1)} & \cdots & \alpha^{(N-1)(N-1)} \end{bmatrix} \begin{bmatrix} x_{0} \\ x_{1} \\ x_{2} \\ \vdots \\ x_{N-1} \end{bmatrix} \]

where \(\alpha = \mathrm{e}^{-\mathrm{i}\frac{2\pi}{N}}\).

Fast Fourier Transform

FFT is an efficient algorithm for implementing DFT, called the Fast Fourier Transform (FFT). It does not make new discoveries in Fourier transform theory, but it greatly advances the application of the Discrete Fourier Transform in computer systems or digital systems. The Number Theoretic Transform (NTT) is an implementation of the Fast Fourier Transform (FFT) based on number theory.

In 1965, Cooley and Tukey published the Fast Fourier Transform algorithm. In fact, FFT had been discovered earlier, but at that time modern computers had not yet emerged, and people did not realize the importance of FFT. Some investigators believe that FFT was discovered by Runge and König in 1924. However, Gauss actually invented this algorithm as early as 1805 but never published it.

Divide and Conquer Implementation

The basic idea of the FFT algorithm is divide and conquer. For DFT, it divides the problem to compute the value of \(f(x)\) when \(x=\omega_n^k\). The divide-and-conquer idea of radix-2 FFT is reflected in splitting the polynomial into odd and even terms.

For example, for a polynomial with \(8\) terms:

\[ f(x) = a_0 + a_1x + a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7 \]

Divide it into two groups according to the parity of the exponents, and then factor out an \(x\) from the right side:

\[ \begin{aligned} f(x) &= (a_0+a_2x^2+a_4x^4+a_6x^6) + (a_1x+a_3x^3+a_5x^5+a_7x^7)\\ &= (a_0+a_2x^2+a_4x^4+a_6x^6) + x(a_1+a_3x^2+a_5x^4+a_7x^6) \end{aligned} \]

Define new functions using the odd and even terms:

\[ \begin{aligned} G(x) &= a_0+a_2x+a_4x^2+a_6x^3\\ H(x) &= a_1+a_3x+a_5x^2+a_7x^3 \end{aligned} \]

Then the original \(f(x)\) can be expressed in terms of the new functions as:

\[ f(x)=G\left(x^2\right) + x \times H\left(x^2\right) \]

Using the property of even-order unit roots \(\omega^i_n = -\omega^{i + n/2}_n\), and the fact that \(G\left(x^2\right)\) and \(H\left(x^2\right)\) are even functions, we know that on the complex plane, the values of \(G(x^2)\) and \(H(x^2)\) at \(\omega^i_n\) and \(\omega^{i+n/2}_n\) are the same. We obtain:

\[ \begin{aligned} f(\omega_n^k) &= G((\omega_n^k)^2) + \omega_n^k \times H((\omega_n^k)^2) \\ &= G(\omega_n^{2k}) + \omega_n^k \times H(\omega_n^{2k}) \\ &= G(\omega_{n/2}^k) + \omega_n^k \times H(\omega_{n/2}^k) \end{aligned} \]

and:

\[ \begin{aligned} f(\omega_n^{k+n/2}) &= G(\omega_n^{2k+n}) + \omega_n^{k+n/2} \times H(\omega_n^{2k+n}) \\ &= G(\omega_n^{2k}) - \omega_n^k \times H(\omega_n^{2k}) \\ &= G(\omega_{n/2}^k) - \omega_n^k \times H(\omega_{n/2}^k) \end{aligned} \]

Therefore, after computing \(G(\omega_{n/2}^k)\) and \(H(\omega_{n/2}^k)\), we can simultaneously obtain \(f(\omega_n^k)\) and \(f(\omega_n^{k+n/2})\). Then we recursively apply DFT to \(G\) and \(H\).

Note that the polynomial length that the divide-and-conquer DFT can handle must be \(2^m(m \in \mathbf{N}^ \ast )\); otherwise, during the division, the left and right sides will have different lengths, and the coefficients on the right side cannot be obtained. So before the first DFT, we need to pad the sequence to a length of \(2^m(m \in \mathbf{N}^\ast )\) (with high-order coefficients padded with \(0\)), making the highest-degree term \(2^m-1\).

When substituting values, because we need to substitute \(n\) different values, we substitute \(\omega_n^0,\omega_n^1,\omega_n^2,\cdots, \omega_n^{n-1} (n=2^m(m \in \mathbf{N}^ \ast ))\), a total of \(2^m\) different values.

Regarding code implementation, the STL provides a template for complex numbers, but it can also be implemented manually. The difference is that using STL's complex allows calling the exp function to compute \(\omega_n\). However, using Euler's formula to compute \(\omega_n\) with imaginary numbers is equivalent.

The above is the introduction to DFT in the FFT algorithm, which transforms a polynomial from coefficient representation to point-value representation.

It is worth noting that because we are using unit complex roots, we need to pad the high-order terms of the \(n\)-term polynomial with zeros so that \(n = 2 ^ k, k \in \mathbf{N}^ \ast\).

Recursive FFT

#include <cmath>
#include <complex>

using Comp = std::complex<double>;  // STL complex

constexpr Comp I(0, 1);  // i
constexpr int MAX_N = 1 << 20;

Comp tmp[MAX_N];

// rev=1, DFT; rev=-1, IDFT
// After applying this function, note the handling of normalization coefficients
void DFT(Comp* f, int n, int rev) {
  if (n == 1) return;
  for (int i = 0; i < n; ++i) tmp[i] = f[i];
  // Even indices on the left, odd indices on the right
  for (int i = 0; i < n; ++i) {
    if (i & 1)
      f[n / 2 + i / 2] = tmp[i];
    else
      f[i / 2] = tmp[i];
  }
  Comp *g = f, *h = f + n / 2;
  // Recursive DFT
  DFT(g, n / 2, rev), DFT(h, n / 2, rev);
  // cur is the current unit root; for k = 0, it corresponds to omega^0_n = 1.
  // step is the difference between two unit roots, i.e., omega^k_n = step*omega^{k-1}*n,
  // defined equivalently as exp(I*(-2*M_PI/n*rev))
  Comp cur(1, 0), step(cos(2 * M_PI / n), sin(-2 * M_PI * rev / n));
  for (int k = 0; k < n / 2;
       ++k) {  // F(omega^k_n) = G(omega^k*{n/2}) + omega^k*n\*H(omega^k*{n/2})
    tmp[k] = g[k] + cur * h[k];
    // F(omega^{k+n/2}*n) = G(omega^k*{n/2}) - omega^k_n*H(omega^k\_{n/2})
    tmp[k + n / 2] = g[k] - cur * h[k];
    cur *= step;
  }
  for (int i = 0; i < n; ++i) f[i] = tmp[i];
}

Time complexity \(O(n\log n)\).

Iterative Implementation

This algorithm can be further optimized from the perspective of "divide and conquer". For radix-2 FFT, we always separate the odd and even coefficients of the entire polynomial until only one coefficient remains. However, this recursive process requires more memory. Therefore, we can first "simulate recursion" to "split" these coefficients in the original array, and then "merge" the computed values in a "doubling" manner.

For "splitting", bit-reversal permutation can be used.

For "merging", butterfly operations can be used to optimize, requiring only \(O(1)\) extra space.

Bit-Reversal Permutation

Take an \(8\)-term polynomial as an example to simulate the splitting process:

• Initial sequence: \(\{x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7\}\)
• After one split: \(\{x_0, x_2, x_4, x_6\},\{x_1, x_3, x_5, x_7 \}\)
• After two splits: \(\{x_0,x_4\} \{x_2, x_6\},\{x_1, x_5\},\{x_3, x_7 \}\)
• After three splits: \(\{x_0\}\{x_4\}\{x_2\}\{x_6\}\{x_1\}\{x_5\}\{x_3\}\{x_7 \}\)

Pattern: Essentially, it is the original sequence, with each number represented in binary, and then the binary representation is reversed symmetrically to obtain the final index. For example, \(x_1\) is 001, reversed to 100, which is 4, and indeed the final position is 4. We call this transformation the bit-reversal permutation; the proof is left to the reader.

According to its definition, we can compute the transformed result for each number in \(O(n)\) time:

Bit-Reversal Permutation Implementation (\(O(n)\))

/*
 * Pre-FFT and IFFT bit-reversal permutation
 * Swap position i with the position of the binary reversal of i
 * len must be a power of 2
 */
void change(Complex y[], int len) {
  // Initially i is 0...01, and j is 10...0, symmetric in binary.
  // Then i increments gradually, while j remains symmetric to i, until i = 1...11.
  for (int i = 1, j = len / 2, k; i < len - 1; i++) {
    // Swap elements with reversed indices; i < j ensures swapping only once
    if (i < j) swap(y[i], y[j]);
    // i does normal +1, j does reversal-type +1, always keeping i and j reversed.
    // Here k represents the highest bit where 0 appears. j first subtracts the high bits all 1,
    // until encountering 0, then adds back.
    // Considering the number of bit flips in j, the highest bit will flip n times,
    // the second highest will flip n/2 times, and so on, so the time complexity is:
    // T(n) = n + n/2 + n/4 + ... = O(n)
    k = len / 2;
    while (j >= k) {
      j = j - k;
      k = k / 2;
    }
    j += k;
  }
}

Bit-reversal permutation can also be implemented iteratively from small to large in \(O(n)\). Let \(len=2^k\), where \(k\) denotes the length of the binary number, and let \(R(x)\) denote the number obtained by reversing the binary digits of \(x\) of length \(k\) (padding with zeros on the high bits). We want to compute \(R(0),R(1),\cdots,R(n-1)\).

First, \(R(0)=0\).

We compute \(R(x)\) from small to large. Therefore, when computing \(R(x)\), the value of \(R\left(\left\lfloor \dfrac{x}{2} \right\rfloor\right)\) is known. So we shift \(x\) right by one bit (divide by \(2\)), then reverse, then shift right by one bit, obtaining the reversal of all bits of \(x\) except the least significant bit.

Consider the reversal of the least significant bit: if the least significant bit is \(0\), the highest bit after reversal is \(0\). If the least significant bit is \(1\), then the highest bit after reversal is \(1\), so we need to add \(\dfrac{len}{2}=2^{k-1}\). In summary:

\[ R(x)=\left\lfloor \frac{R\left(\left\lfloor \frac{x}{2} \right\rfloor\right)}{2} \right\rfloor + (x\bmod 2)\times \frac{len}{2} \]

For example: let \(k=5\), \(len=(100000)_2\). To reverse \((11001)_2\):

1. Consider \((1100)_2\), we know \(R((1100)_2)=R((01100)_2)=(00110)_2\), then shift right by one bit to get \((00011)_2\).
2. Consider the least significant bit: if it is \(1\), it flips to the highest bit of the number, i.e., add \((10000)_2=2^{k-1}\) to the reversed number; if it is \(0\), no change.

Bit-Reversal Permutation Implementation (\(O(n)\))

// Also ensure len is a power of 2
// Let rev[i] be the reversed value of i
void change(Complex y[], int len) {
  for (int i = 0; i < len; ++i) {
    rev[i] = rev[i >> 1] >> 1;
    if (i & 1) {  // If the last bit is 1, flip to len/2
      rev[i] |= len >> 1;
    }
  }
  for (int i = 0; i < len; ++i) {
    if (i < rev[i]) {  // Ensure each pair is swapped only once
      swap(y[i], y[rev[i]]);
    }
  }
  return;
}

Butterfly Operation Optimization

Given \(G(\omega_{n/2}^k)\) and \(H(\omega_{n/2}^k)\), we need to compute \(f(\omega_n^k)\) and \(f(\omega_n^{k+n/2})\) using the following two formulas:

\[ \begin{aligned} f(\omega_n^k) & = G(\omega_{n/2}^k) + \omega_n^k \times H(\omega_{n/2}^k) \\ f(\omega_n^{k+n/2}) & = G(\omega_{n/2}^k) - \omega_n^k \times H(\omega_{n/2}^k) \end{aligned} \]

After applying bit-reversal permutation, for given \(n, k\):

• The value of \(G(\omega_{n/2}^k)\) is stored at array index \(k\), and the value of \(H(\omega_{n/2}^k)\) is stored at array index \(k + \dfrac{n}{2}\).
• The value of \(f(\omega_n^k)\) will be stored at array index \(k\), and the value of \(f(\omega_n^{k+n/2})\) will be stored at array index \(k + \dfrac{n}{2}\).

Therefore, we can directly overwrite at array indices \(k\) and \(k + \frac{n}{2}\) without needing extra arrays to save values. This method is called the butterfly operation, or more precisely, the radix-2 butterfly operation.

To elaborate on how to complete all merge operations with segment length \(\frac{n}{2}\) using butterfly operations:

1. Let segment length \(s = \frac{n}{2}\);
2. Simultaneously enumerate the left endpoints \(l_g = 0, 2s, 4s, \cdots, N-2s\) of the sequence \(\{G(\omega_{n/2}^k)\}\) and the left endpoints \(l_h = s, 3s, 5s, \cdots, N-s\) of the sequence \(\{H(\omega_{n/2}^k)\}\);
3. When merging two segments, enumerate \(k = 0, 1, 2, \cdots, s-1\). At this point, \(G(\omega_{n/2}^k)\) is stored at array index \(l_g + k\), and \(H(\omega_{n/2}^k)\) is stored at array index \(l_h + k\);
4. Use butterfly operations to compute \(f(\omega_n^k)\) and \(f(\omega_n^{k+n/2})\), then directly overwrite in place.

Inverse Fast Fourier Transform

The Inverse Fourier Transform can be expressed using the Fourier Transform. We have two ways to understand this.

Linear Algebra Perspective

The role of IDFT (Inverse Fourier Transform) is to convert the point-value representation of the target polynomial back to coefficient representation. DFT itself is a linear transformation, which can be understood as treating the target polynomial as a vector and left-multiplying by a matrix to obtain the transformed vector, simulating the process of substituting unit complex roots into the polynomial:

\[ \begin{bmatrix}y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{bmatrix} = \begin{bmatrix}1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1} \\ 1 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)} \\ 1 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)^2} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{bmatrix} \]

Now we have obtained the leftmost result. The \(x\) values in the middle correspond one-to-one with the point-value representation of the target polynomial. Therefore, based on basic matrix knowledge, we only need to left-multiply both sides by the inverse of that large matrix.

Since the elements of this matrix are very special, its inverse matrix also has special properties: each element is taken reciprocal and then divided by the transformation length \(n\).

Note: The length of the Fourier Transform is not the length of the polynomial; the transform length should be longer than the length of the product polynomial. If the polynomials to be multiplied are not long enough, pad with \(0\) at high-order terms.

To make the computation result the reciprocal of the original, according to Euler's formula, we have:

\[ \frac{1}{\omega_k}=\omega_k^{-1}=\mathrm{e}^{-\frac{2\pi \mathrm{i}}{k}}=\cos\left(\frac{2\pi}{k}\right)+\mathrm{i} \sin\left(-\frac{2\pi}{k}\right) \]

Therefore, we can try to set the unit root \(\omega_k\) to \(\mathrm{e}^{-\frac{2\pi \mathrm{i}}{k}}\). Then our computation result will become the reciprocal of the original, and the only extra operation is to divide by its length \(n\). All other operations are exactly the same as DFT. We can define a function with an additional parameter \(1\) or \(-1\), and multiply it to \(\pi\). Passing \(1\) means DFT, and passing \(-1\) means IDFT.

Periodicity of Unit Complex Roots

The periodicity of unit complex roots can also be used to understand the relationship between IDFT and DFT.

Consider the original polynomial \(f(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}=\sum_{i=0}^{n-1}a_ix^i\). IDFT converts the point-value representation back to coefficient representation.

Consider a constructive method. We know \(y_i=f\left( \omega_n^i \right),i\in\{0,1,\cdots,n-1\}\), and we want to find \(\{a_0,a_1,\cdots,a_{n-1}\}\). Construct a polynomial as follows:

\[ A(x)=\sum_{i=0}^{n-1}y_ix^i \]

This is equivalent to treating \(\{y_0,y_1,y_2,\cdots,y_{n-1}\}\) as the coefficient representation of polynomial \(A\).

Now we have two derivation methods, corresponding to two implementation approaches.

Method One

Let \(b_i=\omega_n^{-i}\). Then the point-value representation of polynomial \(A\) at \(x=b_0,b_1,\cdots,b_{n-1}\) is \(\left\{ A(b_0),A(b_1),\cdots,A(b_{n-1}) \right\}\).

Transforming the definition of \(A(x)\), we can express \(A(b_k)\) as:

\[ \begin{aligned} A(b_k)&=\sum_{i=0}^{n-1}f(\omega_n^i)\omega_n^{-ik}=\sum_{i=0}^{n-1}\omega_n^{-ik}\sum_{j=0}^{n-1}a_j(\omega_n^i)^{j}\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j\omega_n^{i(j-k)}=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\left(\omega_n^{j-k}\right)^i\\ \end{aligned} \]

Denote \(S\left(\omega_n^a\right)=\sum_{i=0}^{n-1}\left(\omega_n^a\right)^i\).

When \(a=0 \pmod{n}\), \(S\left(\omega_n^a\right)=n\).

When \(a\neq 0 \pmod{n}\), we perform subtraction with offset:

\[ \begin{aligned} S\left(\omega_n^a\right)&=\sum_{i=0}^{n-1}\left(\omega_n^a\right)^i\\ \omega_n^a S\left(\omega_n^a\right)&=\sum_{i=1}^{n}\left(\omega_n^a\right)^i\\ S\left(\omega_n^a\right)&=\frac{\left(\omega_n^a\right)^n-\left(\omega_n^a\right)^0}{\omega_n^a-1}=0\\ \end{aligned} \]

That is:

\[ S\left(\omega_n^a\right)= \begin{cases} n,&a=0\\ 0,&a\neq 0 \end{cases} \]

Substituting back:

\[ A(b_k)=\sum_{j=0}^{n-1}a_jS\left(\omega_n^{j-k}\right)=a_k\cdot n \]

Thus, given points \(b_i=\omega_n^{-i}\), the point-value representation of \(A\) is:

\[ \begin{aligned} &\left\{ (b_0,A(b_0)),(b_1,A(b_1)),\cdots,(b_{n-1},A(b_{n-1})) \right\}\\ =&\left\{ (b_0,a_0\cdot n),(b_1,a_1\cdot n),\cdots,(b_{n-1},a_{n-1}\cdot n) \right\} \end{aligned} \]

In summary, we take the unit roots as their reciprocals, run FFT on \(\{y_0,y_1,y_2,\cdots,y_{n-1}\}\), and then divide by \(n\) to obtain the coefficient representation of \(f(x)\).

Method Two

We directly substitute \(\omega_n^i\) into \(A(x)\).

The derivation process is similar to Method One, and we eventually obtain \(A(\omega_n^k) = \sum_{j=0}^{n-1}a_jS\left(\omega_n^{j+k}\right)\).

If and only if \(j+k=0 \pmod{n}\), \(S\left(\omega_n^{j+k}\right) = n\); otherwise, it is \(0\). Therefore, \(A(\omega_n^k) = a_{n-k}\cdot n\).

This means that after performing DFT on \(\{y_0,y_1,y_2,\cdots,y_{n-1}\}\), dividing by \(n\), and then reversing the last \(n - 1\) elements, we can also recover the coefficient representation of \(f(x)\).

Code Implementation

Thus, our FFT function can combine DFT and IDFT. The code implementation is as follows:

Iterative FFT (corresponding to Method One)

/*
 * Perform FFT
 * len must be of the form 2^k
 * on == 1 for DFT, on == -1 for IDFT
 */
void fft(Complex y[], int len, int on) {
  // Bit-reversal permutation
  change(y, len);
  // Simulate the merging process, starting from length 1 merging to length 2, up to length len.
  for (int h = 2; h <= len; h <<= 1) {
    // wn: current unit root step: w^1_h
    Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
    // Merge, total len / h times.
    for (int j = 0; j < len; j += h) {
      // Compute current unit root, starting from 1 = w^0_n, then increasing with step wn: w^1_n
      // ...
      Complex w(1, 0);
      for (int k = j; k < j + h / 2; k++) {
        // Left part and right part are solutions to subproblems
        Complex u = y[k];
        Complex t = w * y[k + h / 2];
        // This adds the results of the two divide-and-conquer parts
        y[k] = u + t;
        y[k + h / 2] = u - t;
        // The ω in the latter half "step" must be the opposite of that in the former half
        // Points on the "red circle" rotate a full circle "back", rotating half a circle gives the opposite number
        // The square of a number's opposite equals the square of the number itself
        w = w * wn;
      }
    }
  }
  // If IDFT, each element of its inverse matrix is not just the reciprocal of the original element, but also divided by length len.
  if (on == -1) {
    for (int i = 0; i < len; i++) {
      y[i].x /= len;
      y[i].y /= len;
    }
  }
}

Iterative FFT (corresponding to Method Two)

/*
 * Perform FFT
 * len must be of the form 2^k
 * on == 1 for DFT, on == -1 for IDFT
 */
void fft(Complex y[], int len, int on) {
  change(y, len);
  for (int h = 2; h <= len; h <<= 1) {             // Simulate merging process
    Complex wn(cos(2 * PI / h), sin(2 * PI / h));  // Compute current unit root
    for (int j = 0; j < len; j += h) {
      Complex w(1, 0);  // Compute current unit root
      for (int k = j; k < j + h / 2; k++) {
        Complex u = y[k];
        Complex t = w * y[k + h / 2];
        y[k] = u + t;  // This adds the results of the two divide-and-conquer parts
        y[k + h / 2] = u - t;
        // The ω in the latter half "step" must be the opposite of that in the former half
        // Points on the "red circle" rotate a full circle "back", rotating half a circle gives the opposite number
        // The square of a number's opposite equals the square of the number itself
        w = w * wn;
      }
    }
  }
  if (on == -1) {
    reverse(y + 1, y + len);
    for (int i = 0; i < len; i++) {
      y[i].x /= len;
      y[i].y /= len;
    }
  }
}

FFT Template (HDU 1402 - A * B Problem Plus)

#include <cmath>
#include <cstring>
#include <iostream>

const double PI = acos(-1.0);

struct Complex {
  double x, y;

  Complex(double _x = 0.0, double _y = 0.0) {
    x = _x;
    y = _y;
  }

  Complex operator-(const Complex &b) const {
    return Complex(x - b.x, y - b.y);
  }

  Complex operator+(const Complex &b) const {
    return Complex(x + b.x, y + b.y);
  }

  Complex operator*(const Complex &b) const {
    return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
  }
};

/*
 * Bit-reversal permutation before performing FFT and IFFT
 * Swap element at position i with the element at the bit-reversed position of i
 * len must be a power of 2
 */
void change(Complex y[], int len) {
  int i, j, k;

  for (int i = 1, j = len / 2; i < len - 1; i++) {
    if (i < j) std::swap(y[i], y[j]);

    // Swap elements that are bit-reversals of each other; i<j ensures swapping only once
    // i increments normally, j increments in a bit-reversed manner, keeping i and j as bit-reversals
    k = len / 2;

    while (j >= k) {
      j = j - k;
      k = k / 2;
    }

    if (j < k) j += k;
  }
}

/*
 * Perform FFT
 * len must be of the form 2^k
 * on == 1 for DFT, on == -1 for IDFT
 */
void fft(Complex y[], int len, int on) {
  change(y, len);

  for (int h = 2; h <= len; h <<= 1) {
    Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));

    for (int j = 0; j < len; j += h) {
      Complex w(1, 0);

      for (int k = j; k < j + h / 2; k++) {
        Complex u = y[k];
        Complex t = w * y[k + h / 2];
        y[k] = u + t;
        y[k + h / 2] = u - t;
        w = w * wn;
      }
    }
  }

  if (on == -1) {
    for (int i = 0; i < len; i++) {
      y[i].x /= len;
    }
  }
}

constexpr int MAXN = 200020;
Complex x1[MAXN], x2[MAXN];
char str1[MAXN / 2], str2[MAXN / 2];
int sum[MAXN];
using std::cin;
using std::cout;

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  while (cin >> str1 >> str2) {
    int len1 = strlen(str1);
    int len2 = strlen(str2);
    int len = 1;

    while (len < len1 * 2 || len < len2 * 2) len <<= 1;

    for (int i = 0; i < len1; i++) x1[i] = Complex(str1[len1 - 1 - i] - '0', 0);

    for (int i = len1; i < len; i++) x1[i] = Complex(0, 0);

    for (int i = 0; i < len2; i++) x2[i] = Complex(str2[len2 - 1 - i] - '0', 0);

    for (int i = len2; i < len; i++) x2[i] = Complex(0, 0);

    fft(x1, len, 1);
    fft(x2, len, 1);

    for (int i = 0; i < len; i++) x1[i] = x1[i] * x2[i];

    fft(x1, len, -1);

    for (int i = 0; i < len; i++) sum[i] = int(x1[i].x + 0.5);

    for (int i = 0; i < len; i++) {
      sum[i + 1] += sum[i] / 10;
      sum[i] %= 10;
    }

    len = len1 + len2 - 1;

    while (sum[len] == 0 && len > 0) len--;

    for (int i = len; i >= 0; i--) cout << char(sum[i] + '0');

    cout << '\n';
  }

  return 0;
}

References

1. Peach Sauce's Algorithm Notes.

---

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用