Fast Fourier Transform

Prerequisite: Complex numbers.

This article introduces an algorithm that computes the product of two polynomials in $O(n\log n)$ time, where each polynomial has degree $n$, more efficiently than the naive $O(n^2)$ algorithm. Since the multiplication of two integers can also be regarded as polynomial multiplication, this algorithm can accelerate multiplication of large integers.

Introduction

We now introduce two polynomials $A$ and $B$:

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

Clearly, polynomial $C$ has coefficients $c_i$ satisfying $c_i = \sum_{j = 0}^i a_j b_{i - j}$. For this naive algorithm, computing each term takes $O(n)$ time, and there are $O(n)$ terms, so the total time complexity is $O(n^2)$.

Can we accelerate this and reduce its time complexity? Using the fast Fourier transform, we can reduce the complexity to $O(n \log n)$.

Fourier Transform

The Fourier transform is a method for analyzing signals: it can analyze a signal into its components and also synthesize a signal from those components. Many waveforms can serve as signal components; the Fourier transform uses sine waves as its components.

Let $f(t)$ be a function of time $t$. The Fourier transform detects how strongly the period of frequency $\omega$ appears in $f(t)$:

Its inverse transform is

The inverse transform has a form very similar to the forward transform; the denominator $2\pi$ is exactly the period of the exponential function.

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

The Fourier transform has a corresponding convolution theorem, which converts convolution in the time domain into multiplication in the frequency domain, and convolution in the frequency domain into multiplication in the time domain.

Discrete Fourier Transform

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

The Fourier transform is a continuous inner product of functions in integral form, while the discrete Fourier transform is an inner product in summation form.

Let $\{x_n\}_{n=0}^{N-1}$ be a sequence satisfying the finiteness condition. Its discrete Fourier transform (DFT) is:

Here $\mathrm{e}$ is the base of the natural logarithm, and $i$ is the imaginary unit. This transform is usually denoted by $\mathcal {F}$, namely:

Analogously to the integral form, its inverse discrete Fourier transform (IDFT) is:

It can be written as:

In fact, the normalization coefficients before the sums in the DFT and IDFT formulas are not important. In the definition above, the coefficients before the 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. Because of the special nature of the summation form, it admits other interpretations.

If the sequence $x_n$ is regarded as the coefficient of polynomial $f(x)$'s $x^n$ term, then the computed $X_k$ is exactly the value of polynomial $f(x)$ at the root of unity $\mathrm{e}^{\frac{-2\pi \mathrm{i}k}{N}}$, namely $f(\mathrm{e}^{\frac{-2\pi \mathrm{i}k}{N}})$.

This gives another interpretation of the convolution theorem: a special evaluation operation on polynomials. The discrete Fourier transform is exactly the evaluation of a polynomial at roots of unity.

For example, compute:

Define the function $f(x)$ as:

Then we find that substituting the fourth root of unity into $f(\mathrm{i})$ gives the following sequence:

Thus the following sum exactly cancels the other terms:

Therefore, the answer to this mathematical problem is:

This mathematical problem evaluates at roots of unity, which exactly constitutes the discrete Fourier transform.

Matrix formula

Since the discrete Fourier transform is a linear operator, it can be described by matrix multiplication. In matrix form, the discrete Fourier transform is:

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

Fast Fourier Transform

FFT is an efficient algorithm for implementing the DFT, called the fast Fourier transform (Fast Fourier Transform, FFT). It does not make new discoveries about the theory of the Fourier transform, but it was a major step forward for applying the discrete Fourier transform in computer or digital systems. The fast 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 modern computers did not yet exist, so people did not realize its importance. Some researchers believe FFT was discovered by Runge and König in 1924. In fact, Gauss had invented the 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 the DFT, it recursively computes the value at $x=\omega_n^k$ of $f(x)$. The divide-and-conquer idea of radix-2 FFT is to process the odd-degree and even-degree terms of a polynomial separately.

For example, consider a polynomial with $8$ terms:

Split the terms into two groups according to the parity of their degrees, and factor out an $x$ from the right group:

Construct new functions from the odd- and even-degree terms separately:

Then the original $f(x)$ can be represented using the new functions as:

Using the property of roots of unity $\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 at $\omega^i_n$ and $\omega^{i+n/2}_n$, the corresponding values of $G(x^2)$ and $H(x^2)$ on the complex plane are the same. We obtain:

and:

Thus, after computing $G(\omega_{n/2}^k)$ and $H(\omega_{n/2}^k)$, we can simultaneously compute $f(\omega_n^k)$ and $f(\omega_n^{k+n/2})$. We then recursively perform the DFT on $G$ and $H$, respectively.

Since divide-and-conquer DFT can only handle polynomials whose lengths are $2^m(m \in \mathbf{N}^ \ast )$, otherwise the two sides have different lengths during the division and the right side has no coefficient to access. Therefore, before the first DFT, pad the sequence upward to length $2^m(m \in \mathbf{N}^\ast )$ (padding high-degree coefficients with $0$), so that the highest degree is $2^m-1$.

When substituting values, because we need to substitute $n$ distinct 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$ distinct values.

For implementation, the STL provides a complex-number template, though it can also be implemented manually. The difference is that STL complex can call the exp function to obtain $\omega_n$. In fact, using the imaginary number obtained from Euler's formula to obtain $\omega_n$ is equivalent.

This completes the introduction to the DFT part of the FFT algorithm: it converts a polynomial from coefficient representation to point-value representation.

It is worth noting that because these are roots of unity, we need to pad the high-order part of an $n$-term polynomial with zeros so that $n = 2 ^ k, k \in \mathbf{N}^ \ast$.

#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
// Pay attention to the normalization coefficient after applying this function
void DFT(Comp* f, int n, int rev) {
  if (n == 1) return;
  for (int i = 0; i < n; ++i) tmp[i] = f[i];
  // Put even terms on the left and odd terms 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 root of unity; for k = 0, the corresponding root of unity is omega^0_n = 1.
  // step is the ratio between two consecutive roots of unity, satisfying omega^k_n = step*omega^{k-1}*n,
  // which is equivalently defined 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];
}

The time complexity is $O(n\log n)$.

Iterative doubling implementation

This algorithm can be further optimized from the perspective of divide and conquer. For radix-2 FFT, at each step we separate the odd- and even-degree coefficients of the entire polynomial, continuing until only one coefficient remains. However, this recursive process requires more memory. Therefore, we can first “split” these coefficients in the original array by “mimicking recursion”, and then “double” the segment length to merge the computed values.

For the “split”, we can use a bit-reversal permutation.

For the “merge”, butterfly-operation optimization can complete the work using only $O(1)$ extra space.

Bit-reversal permutation

Using an $8$-term polynomial as an example, simulate the splitting process:

• The initial sequence is $\{x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7\}$.\n- After one bisection: $\{x_0, x_2, x_4, x_6\},\{x_1, x_3, x_5, x_7 \}$.\n- After two bisections: $\{x_0,x_4\} \{x_2, x_6\},\{x_1, x_5\},\{x_3, x_7 \}$.\n- After three bisections: $\{x_0\}\{x_4\}\{x_2\}\{x_6\}\{x_1\}\{x_5\}\{x_3\}\{x_7 \}$.

The pattern is as follows: in the original sequence, write each number in binary and reverse the binary digits symmetrically; the result is the index of its final position. For example, $x_1$ is 001, whose reversal is 100, namely 4, and its final position is indeed 4. We call this transformation the bit-reversal permutation; the proof is left to the reader.

By its definition, we can compute the transformed position of every number in $O(n)$ time:

/*
 * Perform the reversal transformation before FFT and IFFT
 * Swap position i with the position obtained by reversing the binary digits 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; they are symmetric opposites in binary.
  // Then i gradually increases by one while j remains the binary reverse of i, until i = 1...11.
  for (int i = 1, j = len / 2, k; i < len - 1; i++) {
    // Swap elements whose indices are reverses of each other; i < j ensures each pair is swapped once
    if (i < j) swap(y[i], y[j]);
    // i performs an ordinary +1, while j performs a reversal-type +1, keeping i and j as reverses.
    // Here k denotes the highest position at which 0 appears. First subtract the number whose high bits are all 1 from j,
    // until reaching a 0, then add it back.
    // Considering the number of bit flips in j, the highest bit flips n times,
    // the second-highest bit flips 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;
  }
}

The bit-reversal permutation can also be computed recursively in increasing order in $O(n)$ time. Let $len=2^k$, where $k$ is the length of the binary numbers, and let $R(x)$ denote the number obtained by reversing the $k$-digit binary representation of $x$ (padding with leading $0$s). We need $R(0),R(1),\cdots,R(n-1)$.

First, $R(0)=0$.

We compute $R(x)$ in increasing order. Thus, when computing $R(x)$, the value of $R\left(\left\lfloor \dfrac{x}{2} \right\rfloor\right)$ is known. We shift $x$ right by one bit (divide by $2$), then reverse it and shift right by one bit again; this gives the reversal of all bits of $x$ except the least significant (binary) bit.

Consider the reversal of the least significant bit: if it is $0$, the highest bit after reversal is $0$. If it is $1$, the highest bit after reversal is $1$, so we must add $\dfrac{len}{2}=2^{k-1}$. Therefore:

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

1. Consider $(1100)_2$. We know $R((1100)_2)=R((01100)_2)=(00110)_2$; shifting right by one bit gives $(00011)_2$.\n2. Consider the least significant bit. If it is $1$, it is reversed to the highest bit, so add $(10000)_2=2^{k-1}$ to the reversed number; if it is $0$, make no change.

// Again, ensure that len is a power of 2
// Let rev[i] be the value obtained by reversing 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, reverse it to len/2
      rev[i] |= len >> 1;
    }
  }
  for (int i = 0; i < len; ++i) {
    if (i < rev[i]) {  // ensure that each pair is reversed only once
      swap(y[i], y[rev[i]]);
    }
  }
  return;
}

Butterfly-operation optimization

After obtaining $G(\omega_{n/2}^k)$ and $H(\omega_{n/2}^k)$, use the following two formulas to compute $f(\omega_n^k)$ and $f(\omega_n^{k+n/2})$:

After applying the bit-reversal permutation, for a given pair $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}$.\n- 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 overwrite the values directly at array indices $k$ and $k + \frac{n}{2}$ without allocating an extra array to save them. This method is called the butterfly operation, or more precisely, the radix-2 butterfly operation.

Here is a more detailed explanation of how to use butterfly operations to merge all segments of length $\frac{n}{2}$:

1. Let the segment length be $s = \frac{n}{2}$.\n2. Enumerate simultaneously the sequence $\{G(\omega_{n/2}^k)\}$ and its left endpoint $l_g = 0, 2s, 4s, \cdots, N-2s$, and the sequence $\{H(\omega_{n/2}^k)\}$ and its left endpoint $l_h = s, 3s, 5s, \cdots, N-s$.\n3. When merging two segments, enumerate $k = 0, 1, 2, \cdots, s-1$; at this time, $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$.\n4. Use the butterfly operation to compute $f(\omega_n^k)$ and $f(\omega_n^{k+n/2})$, then overwrite the original positions directly.

Inverse Fast Fourier Transform

The inverse Fourier transform can be expressed using the Fourier transform. There are two ways to understand this.

From the perspective of linear algebra

The IDFT (inverse Fourier transform) converts the point-value representation of the target polynomial into its coefficient representation. Since the DFT itself is a linear transformation, we can regard the target polynomial as a vector and left-multiply it by a matrix to obtain the transformed vector, simulating the process of substituting roots of unity into the polynomial:

We have now obtained the result on the far left. The middle $x$ values also correspond one-to-one in the point-value representation of the target polynomial. Therefore, by basic matrix theory, we only need to left-multiply both sides of the equation by the inverse of the large matrix in the middle.

Because the entries of this matrix are very special, its inverse also has a special property: taking the reciprocal of each entry and then dividing by the transform length $n$ gives its inverse matrix.

Note that 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 their high-degree terms with $0$.

To make the computed result the reciprocal of the original, Euler's formula gives:

Therefore, we can try taking the root of unity $\omega_k$ to be $\mathrm{e}^{-\frac{2\pi \mathrm{i}}{k}}$. This makes the computed result the reciprocal of the original; the only extra operation is then to divide by its length $n$, while all other operations are exactly the same as in the DFT. We can define a function with a parameter $1$ or $-1$ multiplied by $\pi$. Passing $1$ means DFT, and passing $-1$ means IDFT.

Periodicity of roots of unity

The periodicity of roots of unity gives another way 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$. The IDFT restores the coefficient representation from the point-value representation.

Consider the construction method. We know $y_i=f\left( \omega_n^i \right),i\in\{0,1,\cdots,n-1\}$ and want to find $\{a_0,a_1,\cdots,a_{n-1}\}$. Construct the following polynomial:

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

At this point, there are two derivations, corresponding to two implementation methods.

Method 1

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\}$.

By transforming the definition of $A(x)$, we can express $A(b_k)$ as

Let $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}$, subtract the equations after shifting one term:

That is,

Substituting back into the original equation gives

In other words, given the points $b_i=\omega_n^{-i}$, the point-value representation of $A$ is

In summary, take the reciprocal of the root of unity, run FFT on $\{y_0,y_1,y_2,\cdots,y_{n-1}\}$, and divide by $n$ to obtain the coefficient representation of $f(x)$.

Method 2

Directly substitute $\omega_n^i$ into $A(x)$.

The derivation is much the same as in Method 1; ultimately we obtain $A(\omega_n^k) = \sum_{j=0}^{n-1}a_jS\left(\omega_n^{j+k}\right)$.

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

This means that after applying the DFT to $\{y_0,y_1,y_2,\cdots,y_{n-1}\}$ and dividing by $n$, by reversing the last $n - 1$ elements we can restore the coefficient representation of $f(x)$.

Code implementation

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

/*
 * Perform FFT
 * len must have the form 2^k
 * on == 1 means DFT, and on == -1 means IDFT
 */
void fft(Complex y[], int len, int on) {
  // Bit-reversal permutation
  change(y, len);
  // Simulate the merge process, starting by merging length one into length two and continuing up to length len.
  for (int h = 2; h <= len; h <<= 1) {
    // wn: the interval between the current roots of unity: w^1_h
    Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
    // Merge, for a total of len / h times.
    for (int j = 0; j < len; j += h) {
      // Compute the current root of unity; initially it is 1 = w^0_n, then it increases with wn as the interval: w^1_n
      // ...
      Complex w(1, 0);
      for (int k = j; k < j + h / 2; k++) {
        // The left and right parts are the solutions to the 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 second half of the “step” is necessarily the opposite of the one in the “first half”
        // A point on the “red circle” returns after one full turn and becomes its opposite after a half turn
        // The square of a number's opposite equals the square of the number itself
        w = w * wn;
      }
    }
  }
  // For IDFT, each element of its inverse matrix is not only the reciprocal of the original element; it must also be divided by the length len.
  if (on == -1) {
    for (int i = 0; i < len; i++) {
      y[i].x /= len;
      y[i].y /= len;
    }
  }
}

/*
 * Perform FFT
 * len must have the form 2^k
 * on == 1 means DFT, and on == -1 means IDFT
 */
void fft(Complex y[], int len, int on) {
  change(y, len);
  for (int h = 2; h <= len; h <<= 1) {             // Simulate the merge process
    Complex wn(cos(2 * PI / h), sin(2 * PI / h));  // Compute the current root of unity
    for (int j = 0; j < len; j += h) {
      Complex w(1, 0);  // Compute the current root of unity
      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 second half of the “step” is necessarily the opposite of the one in the “first half”
        // A point on the “red circle” returns after one full turn and becomes its opposite after a half turn
        // 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;
    }
  }
}

#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);
  }
};

/*
 * Perform the reversal transformation before FFT and IFFT
 * Swap position i with the position obtained by reversing the binary digits 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 whose indices are reverses of each other; i<j ensures each pair is swapped once
    // i performs an ordinary +1, while j performs a reversal-type +1, keeping i and j as reverses
    k = len / 2;

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

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

/*
 * Perform FFT
 *len must have the form 2^k
 *on == 1 means DFT, and on == -1 means 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. Tao Jiang's Algorithm Notes.