cp-library

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub weilycoder/cp-library

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

Depends on

Verified with

Code

#ifndef WEILYCODER_POLY_FFT_CONVOLVE_HPP
#define WEILYCODER_POLY_FFT_CONVOLVE_HPP

#include "fft.hpp"
#include <complex>
#include <cstddef>
#include <vector>

/**
 * @file fft_convolve.hpp
 * @brief Functions for performing convolution using Fast Fourier Transform (FFT).  
 */

namespace weilycoder {
/**
 * @brief Perform convolution of two complex vectors using FFT.
 * @tparam float_t Floating-point type for computations (default: double).
 * @param a First input vector of complex numbers.
 * @param b Second input vector of complex numbers.
 * @return Vector containing the convolution result.
 */
template <typename float_t = double>
std::vector<std::complex<float_t>> fft_convolve(std::vector<std::complex<float_t>> a,
                                                std::vector<std::complex<float_t>> b) {
  size_t n = 1;
  while (n < a.size() + b.size() - 1)
    n <<= 1;
  a.resize(n), b.resize(n);
  fft(a), fft(b);
  for (size_t i = 0; i < n; ++i)
    a[i] *= b[i];
  fft<-1>(a), a.resize(a.size() + b.size() - 1);
  return a;
}

/**
 * @brief Perform convolution of two real-valued vectors using FFT.
 * @tparam float_t Floating-point type for computations (default: double).
 * @param a First input vector of real numbers.
 * @param b Second input vector of real numbers.
 * @return Vector containing the convolution result.
 * @note This function uses a technique that combines two real sequences into
 *       a single complex sequence to perform the convolution more efficiently.
 */
template <typename float_t = double>
std::vector<float_t> fft_convolve_real(const std::vector<float_t> &a,
                                       const std::vector<float_t> &b) {
  size_t n = 1;
  while (n < a.size() + b.size() - 1)
    n <<= 1;
  std::vector<std::complex<float_t>> F(n);
  for (size_t i = 0; i < a.size(); ++i)
    F[i].real(a[i]);
  for (size_t i = 0; i < b.size(); ++i)
    F[i].imag(b[i]);
  fft(F);
  for (size_t i = 0; i < n; ++i)
    F[i] *= F[i];
  fft<-1>(F);
  std::vector<float_t> result(a.size() + b.size() - 1);
  for (size_t i = 0; i < result.size(); ++i)
    result[i] = F[i].imag() / 2;
  return result;
}
} // namespace weilycoder

#endif
#line 1 "weilycoder/poly/fft_convolve.hpp"



#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


#line 8 "weilycoder/poly/fft_convolve.hpp"

/**
 * @file fft_convolve.hpp
 * @brief Functions for performing convolution using Fast Fourier Transform (FFT).  
 */

namespace weilycoder {
/**
 * @brief Perform convolution of two complex vectors using FFT.
 * @tparam float_t Floating-point type for computations (default: double).
 * @param a First input vector of complex numbers.
 * @param b Second input vector of complex numbers.
 * @return Vector containing the convolution result.
 */
template <typename float_t = double>
std::vector<std::complex<float_t>> fft_convolve(std::vector<std::complex<float_t>> a,
                                                std::vector<std::complex<float_t>> b) {
  size_t n = 1;
  while (n < a.size() + b.size() - 1)
    n <<= 1;
  a.resize(n), b.resize(n);
  fft(a), fft(b);
  for (size_t i = 0; i < n; ++i)
    a[i] *= b[i];
  fft<-1>(a), a.resize(a.size() + b.size() - 1);
  return a;
}

/**
 * @brief Perform convolution of two real-valued vectors using FFT.
 * @tparam float_t Floating-point type for computations (default: double).
 * @param a First input vector of real numbers.
 * @param b Second input vector of real numbers.
 * @return Vector containing the convolution result.
 * @note This function uses a technique that combines two real sequences into
 *       a single complex sequence to perform the convolution more efficiently.
 */
template <typename float_t = double>
std::vector<float_t> fft_convolve_real(const std::vector<float_t> &a,
                                       const std::vector<float_t> &b) {
  size_t n = 1;
  while (n < a.size() + b.size() - 1)
    n <<= 1;
  std::vector<std::complex<float_t>> F(n);
  for (size_t i = 0; i < a.size(); ++i)
    F[i].real(a[i]);
  for (size_t i = 0; i < b.size(); ++i)
    F[i].imag(b[i]);
  fft(F);
  for (size_t i = 0; i < n; ++i)
    F[i] *= F[i];
  fft<-1>(F);
  std::vector<float_t> result(a.size() + b.size() - 1);
  for (size_t i = 0; i < result.size(); ++i)
    result[i] = F[i].imag() / 2;
  return result;
}
} // namespace weilycoder
Back to top page