Floating-Point Convolution
(math/fps/floating_point_convolution.hpp)

• View this file on GitHub
• Last update: 2026-07-07 14:26:59+09:00
• Include: #include "math/fps/floating_point_convolution.hpp"

Overview

This header computes polynomial convolution using the complex fast Fourier transform.

It supports real coefficients of type float, double, or long double, and complex coefficients using std::complex of those types.

Functions

convolution_fft(first, second);

For real inputs, the result is a vector of the same real type. For complex inputs, the result is a vector of the same complex type.

convolution_round(first, second);

For integral inputs, this performs a long double FFT and rounds every coefficient to the nearest integer of the input type.

If either input is empty, the result is empty. Otherwise, its size is first.size() + second.size() - 1.

Complexity

Let $N$ be the output size rounded up to a power of two.

• Time: $O(N\log N)$
• Additional memory: $O(N)$

Accuracy

FFT convolution is approximate. Error grows with transform length and coefficient magnitude.

convolution_round is correct only when numerical error remains below one half and every exact coefficient is representable by the output integer type. For exact modular convolution, use math/fps/convolution.hpp. For exact signed long long convolution, use math/fps/convolution_ll.hpp.

Using long double generally improves accuracy, although its precision is platform-dependent.

Example

#include "math/fps/floating_point_convolution.hpp"

#include <iostream>
#include <vector>

int main() {
    std::vector<double> first = {1.5, 2.0};
    std::vector<double> second = {3.0, -1.0};

    auto result = m1une::fps::convolution_fft(first, second);
    for (double value : result) std::cout << value << "\n";
}

Required by

• Math All (math/all.hpp)
• Formal Power Series All (math/fps/all.hpp)

Verified with

• verify/math/fps/floating_point_convolution.test.cpp
• verify/math/fps/fps_algorithms.test.cpp
• verify/math/math_algorithms.test.cpp

Code

#ifndef M1UNE_FPS_FLOATING_POINT_CONVOLUTION_HPP
#define M1UNE_FPS_FLOATING_POINT_CONVOLUTION_HPP 1

#include <algorithm>
#include <bit>
#include <cmath>
#include <complex>
#include <concepts>
#include <numbers>
#include <type_traits>
#include <vector>

namespace m1une {
namespace fps {

namespace floating_point_convolution_detail {

template <std::floating_point Real>
void fft(std::vector<std::complex<Real>>& values, bool inverse) {
    int size = int(values.size());
    for (int index = 1, reversed = 0; index < size; ++index) {
        int bit = size >> 1;
        while (reversed & bit) {
            reversed ^= bit;
            bit >>= 1;
        }
        reversed ^= bit;
        if (index < reversed) std::swap(values[index], values[reversed]);
    }

    for (int length = 2; length <= size; length <<= 1) {
        Real angle = Real(2) * std::numbers::pi_v<Real> / Real(length);
        if (inverse) angle = -angle;
        std::complex<Real> step(std::cos(angle), std::sin(angle));
        int half = length >> 1;
        for (int offset = 0; offset < size; offset += length) {
            std::complex<Real> root(1, 0);
            for (int index = 0; index < half; ++index) {
                std::complex<Real> even = values[offset + index];
                std::complex<Real> odd = values[offset + index + half] * root;
                values[offset + index] = even + odd;
                values[offset + index + half] = even - odd;
                root *= step;
            }
        }
    }

    if (inverse) {
        Real inverse_size = Real(1) / Real(size);
        for (auto& value : values) value *= inverse_size;
    }
}

template <std::floating_point Real>
std::vector<std::complex<Real>> complex_convolution(const std::vector<std::complex<Real>>& first,
                                                    const std::vector<std::complex<Real>>& second) {
    if (first.empty() || second.empty()) return {};
    std::size_t result_size = first.size() + second.size() - 1;
    std::size_t size = std::bit_ceil(result_size);
    std::vector<std::complex<Real>> transformed_first(size);
    std::vector<std::complex<Real>> transformed_second(size);
    std::copy(first.begin(), first.end(), transformed_first.begin());
    std::copy(second.begin(), second.end(), transformed_second.begin());

    fft(transformed_first, false);
    fft(transformed_second, false);
    for (std::size_t index = 0; index < size; ++index) {
        transformed_first[index] *= transformed_second[index];
    }
    fft(transformed_first, true);
    transformed_first.resize(result_size);
    return transformed_first;
}

}  // namespace floating_point_convolution_detail

// Convolution of complex floating-point coefficients.
template <std::floating_point Real>
std::vector<std::complex<Real>> convolution_fft(const std::vector<std::complex<Real>>& first,
                                                const std::vector<std::complex<Real>>& second) {
    return floating_point_convolution_detail::complex_convolution(first, second);
}

// Convolution of real floating-point coefficients.
template <std::floating_point Real>
std::vector<Real> convolution_fft(const std::vector<Real>& first, const std::vector<Real>& second) {
    if (first.empty() || second.empty()) return {};
    std::vector<std::complex<Real>> complex_first(first.size());
    std::vector<std::complex<Real>> complex_second(second.size());
    for (std::size_t index = 0; index < first.size(); ++index) {
        complex_first[index] = std::complex<Real>(first[index], 0);
    }
    for (std::size_t index = 0; index < second.size(); ++index) {
        complex_second[index] = std::complex<Real>(second[index], 0);
    }
    auto result = floating_point_convolution_detail::complex_convolution(complex_first, complex_second);
    std::vector<Real> real_result(result.size());
    for (std::size_t index = 0; index < result.size(); ++index) {
        real_result[index] = result[index].real();
    }
    return real_result;
}

// Uses long-double FFT and rounds each coefficient to the nearest integer.
template <std::integral Integer>
std::vector<Integer> convolution_round(const std::vector<Integer>& first, const std::vector<Integer>& second) {
    if (first.empty() || second.empty()) return {};
    std::vector<long double> real_first(first.begin(), first.end());
    std::vector<long double> real_second(second.begin(), second.end());
    std::vector<long double> real_result = convolution_fft(real_first, real_second);
    std::vector<Integer> result(real_result.size());
    for (std::size_t index = 0; index < result.size(); ++index) {
        result[index] = static_cast<Integer>(std::round(real_result[index]));
    }
    return result;
}

}  // namespace fps
}  // namespace m1une

#endif  // M1UNE_FPS_FLOATING_POINT_CONVOLUTION_HPP

#line 1 "math/fps/floating_point_convolution.hpp"



#include <algorithm>
#include <bit>
#include <cmath>
#include <complex>
#include <concepts>
#include <numbers>
#include <type_traits>
#include <vector>

namespace m1une {
namespace fps {

namespace floating_point_convolution_detail {

template <std::floating_point Real>
void fft(std::vector<std::complex<Real>>& values, bool inverse) {
    int size = int(values.size());
    for (int index = 1, reversed = 0; index < size; ++index) {
        int bit = size >> 1;
        while (reversed & bit) {
            reversed ^= bit;
            bit >>= 1;
        }
        reversed ^= bit;
        if (index < reversed) std::swap(values[index], values[reversed]);
    }

    for (int length = 2; length <= size; length <<= 1) {
        Real angle = Real(2) * std::numbers::pi_v<Real> / Real(length);
        if (inverse) angle = -angle;
        std::complex<Real> step(std::cos(angle), std::sin(angle));
        int half = length >> 1;
        for (int offset = 0; offset < size; offset += length) {
            std::complex<Real> root(1, 0);
            for (int index = 0; index < half; ++index) {
                std::complex<Real> even = values[offset + index];
                std::complex<Real> odd = values[offset + index + half] * root;
                values[offset + index] = even + odd;
                values[offset + index + half] = even - odd;
                root *= step;
            }
        }
    }

    if (inverse) {
        Real inverse_size = Real(1) / Real(size);
        for (auto& value : values) value *= inverse_size;
    }
}

template <std::floating_point Real>
std::vector<std::complex<Real>> complex_convolution(const std::vector<std::complex<Real>>& first,
                                                    const std::vector<std::complex<Real>>& second) {
    if (first.empty() || second.empty()) return {};
    std::size_t result_size = first.size() + second.size() - 1;
    std::size_t size = std::bit_ceil(result_size);
    std::vector<std::complex<Real>> transformed_first(size);
    std::vector<std::complex<Real>> transformed_second(size);
    std::copy(first.begin(), first.end(), transformed_first.begin());
    std::copy(second.begin(), second.end(), transformed_second.begin());

    fft(transformed_first, false);
    fft(transformed_second, false);
    for (std::size_t index = 0; index < size; ++index) {
        transformed_first[index] *= transformed_second[index];
    }
    fft(transformed_first, true);
    transformed_first.resize(result_size);
    return transformed_first;
}

}  // namespace floating_point_convolution_detail

// Convolution of complex floating-point coefficients.
template <std::floating_point Real>
std::vector<std::complex<Real>> convolution_fft(const std::vector<std::complex<Real>>& first,
                                                const std::vector<std::complex<Real>>& second) {
    return floating_point_convolution_detail::complex_convolution(first, second);
}

// Convolution of real floating-point coefficients.
template <std::floating_point Real>
std::vector<Real> convolution_fft(const std::vector<Real>& first, const std::vector<Real>& second) {
    if (first.empty() || second.empty()) return {};
    std::vector<std::complex<Real>> complex_first(first.size());
    std::vector<std::complex<Real>> complex_second(second.size());
    for (std::size_t index = 0; index < first.size(); ++index) {
        complex_first[index] = std::complex<Real>(first[index], 0);
    }
    for (std::size_t index = 0; index < second.size(); ++index) {
        complex_second[index] = std::complex<Real>(second[index], 0);
    }
    auto result = floating_point_convolution_detail::complex_convolution(complex_first, complex_second);
    std::vector<Real> real_result(result.size());
    for (std::size_t index = 0; index < result.size(); ++index) {
        real_result[index] = result[index].real();
    }
    return real_result;
}

// Uses long-double FFT and rounds each coefficient to the nearest integer.
template <std::integral Integer>
std::vector<Integer> convolution_round(const std::vector<Integer>& first, const std::vector<Integer>& second) {
    if (first.empty() || second.empty()) return {};
    std::vector<long double> real_first(first.begin(), first.end());
    std::vector<long double> real_second(second.begin(), second.end());
    std::vector<long double> real_result = convolution_fft(real_first, real_second);
    std::vector<Integer> result(real_result.size());
    for (std::size_t index = 0; index < result.size(); ++index) {
        result[index] = static_cast<Integer>(std::round(real_result[index]));
    }
    return result;
}

}  // namespace fps
}  // namespace m1une

Back to top page