Implementation of Fast Fourier Transform (FFT) and its inverse.
(weilycoder/poly/fft.hpp)

• View this file on GitHub
• Last update: 2025-11-07 09:14:02+08:00
• Include: #include "weilycoder/poly/fft.hpp"

Depends on

• Utility functions and constants for Fast Fourier Transform (FFT) (weilycoder/poly/fft_utility.hpp)

Required by

• Functions for performing convolution using Fast Fourier Transform (FFT). (weilycoder/poly/fft_convolve.hpp)

Verified with

• test/multiplication_of_big_integers.fft.test.cpp

Code

#ifndef WEILYCODER_POLY_FFT_HPP
#define WEILYCODER_POLY_FFT_HPP

#include "fft_utility.hpp"
#include <complex>
#include <cstddef>
#include <stdexcept>
#include <vector>

/**
 * @file fft.hpp
 * @brief Implementation of Fast Fourier Transform (FFT) and its inverse.
 */

namespace weilycoder {
/**
 * @brief In-place iterative Cooley–Tukey FFT / inverse FFT on a complex vector.
 *        The length of the input vector should be a power of two.
 * @tparam on Direction of the transform: 1 for FFT, -1 for inverse FFT.
 * @tparam float_t Floating-point type for computations (default: double).
 * @param y Vector of complex numbers to be transformed in-place.
 */
template <int32_t on = 1, typename float_t = double>
void fft(std::vector<std::complex<float_t>> &y) {
  static_assert(on == 1 || on == -1, "on must be 1 or -1");
  fft_change(y);
  size_t len = y.size();
  if (len == 0 || (len & (len - 1)) != 0)
    throw std::invalid_argument("Length of input vector must be a power of two");
  for (size_t h = 2; h <= len; h <<= 1) {
    std::complex<float_t> wn(cos(2 * PI<float_t> / h), sin(on * 2 * PI<float_t> / h));
    for (size_t j = 0; j < len; j += h) {
      std::complex<float_t> w(1, 0);
      for (size_t k = j; k < j + (h >> 1); ++k, w *= wn) {
        std::complex<float_t> u = y[k], t = w * y[k + (h >> 1)];
        y[k] = u + t, y[k + (h >> 1)] = u - t;
      }
    }
  }
  if constexpr (on == -1)
    for (size_t i = 0; i < len; ++i)
      y[i] /= len;
}
} // namespace weilycoder

#endif

#line 1 "weilycoder/poly/fft.hpp"



#line 1 "weilycoder/poly/fft_utility.hpp"



#include <complex>
#include <cstddef>
#include <vector>

/**
 * @file fft_utility.hpp
 * @brief Utility functions and constants for Fast Fourier Transform (FFT)
 */

namespace weilycoder {
/**
 * @brief Alias for the commonly used complex number type with double precision.
 */
using complex_t = std::complex<double>;

/**
 * @brief Compile-time constant for π (pi) parameterized by numeric type.
 * @tparam T Numeric type (e.g., float, double, long double).
 */
template <typename T> constexpr T PI = T(3.1415926535897932384626433832795l);

/**
 * @brief Perform in-place bit-reversal permutation (index reordering) required by
 *        iterative FFT.
 * @tparam T Element type stored in the vector. Must be Swappable and efficiently
 *         copyable.
 * @param a Vector to be permuted in-place. Its size should typically be a power
 *          of two for use with the FFT implementation in this file.
 */
template <typename T> void fft_change(std::vector<T> &a) {
  size_t n = a.size();
  std::vector<size_t> rev(n);
  for (size_t i = 0; i < n; ++i) {
    rev[i] = rev[i >> 1] >> 1;
    if (i & 1)
      rev[i] |= n >> 1;
    if (i < rev[i])
      std::swap(a[i], a[rev[i]]);
  }
}
} // namespace weilycoder


#line 7 "weilycoder/poly/fft.hpp"
#include <stdexcept>
#line 9 "weilycoder/poly/fft.hpp"

/**
 * @file fft.hpp
 * @brief Implementation of Fast Fourier Transform (FFT) and its inverse.
 */

namespace weilycoder {
/**
 * @brief In-place iterative Cooley–Tukey FFT / inverse FFT on a complex vector.
 *        The length of the input vector should be a power of two.
 * @tparam on Direction of the transform: 1 for FFT, -1 for inverse FFT.
 * @tparam float_t Floating-point type for computations (default: double).
 * @param y Vector of complex numbers to be transformed in-place.
 */
template <int32_t on = 1, typename float_t = double>
void fft(std::vector<std::complex<float_t>> &y) {
  static_assert(on == 1 || on == -1, "on must be 1 or -1");
  fft_change(y);
  size_t len = y.size();
  if (len == 0 || (len & (len - 1)) != 0)
    throw std::invalid_argument("Length of input vector must be a power of two");
  for (size_t h = 2; h <= len; h <<= 1) {
    std::complex<float_t> wn(cos(2 * PI<float_t> / h), sin(on * 2 * PI<float_t> / h));
    for (size_t j = 0; j < len; j += h) {
      std::complex<float_t> w(1, 0);
      for (size_t k = j; k < j + (h >> 1); ++k, w *= wn) {
        std::complex<float_t> u = y[k], t = w * y[k + (h >> 1)];
        y[k] = u + t, y[k + (h >> 1)] = u - t;
      }
    }
  }
  if constexpr (on == -1)
    for (size_t i = 0; i < len; ++i)
      y[i] /= len;
}
} // namespace weilycoder

Back to top page