verify/aoj-ITP1_1_A-discrete_fourier_transform.test.cpp

• View this file on GitHub
• Last update: 2025-07-03 00:41:25+09:00
• Problem: https://onlinejudge.u-aizu.ac.jp/courses/lesson/2/ITP1/1/ITP1_1_A

Depends on

• Discrete Fourier Transform（離散フーリエ変換） (algorithm/Math/Convolution/discrete_fourier_transform.hpp)

Code

#define PROBLEM "https://onlinejudge.u-aizu.ac.jp/courses/lesson/2/ITP1/1/ITP1_1_A"

#include <cassert>
#include <cmath>
#include <iostream>
#include <random>
#include <vector>

#include "../algorithm/Math/Convolution/discrete_fourier_transform.hpp"

int main() {
    constexpr int t = 100;
    std::random_device seed;
    std::mt19937_64 rng(seed());

    constexpr double eps = 1e-6;
    std::uniform_int_distribution<int> uniform_n(1, 200);
    std::uniform_real_distribution<double> uniform(-1e3, 1e3);
    for(int i = 0; i < t; ++i) {
        const int n = uniform_n(rng), m = uniform_n(rng);
        std::vector<double> a(n), b(m);
        for(auto &elem : a) elem = uniform(rng);
        for(auto &elem : b) elem = uniform(rng);

        auto &&target = algorithm::dft::convolve(a, b);
        auto &&want = algorithm::dft::convolve_naive(a, b);

        assert(target.size() == size_t(n + m - 1));
        assert(want.size() == size_t(n + m - 1));
        for(int j = 0; j < n + m - 1; ++j) assert(std::abs(target[j] - want[j]) < eps);
    }

    std::cout << "Hello World" << std::endl;
}

#line 1 "verify/aoj-ITP1_1_A-discrete_fourier_transform.test.cpp"
#define PROBLEM "https://onlinejudge.u-aizu.ac.jp/courses/lesson/2/ITP1/1/ITP1_1_A"

#include <cassert>
#include <cmath>
#include <iostream>
#include <random>
#include <vector>

#line 1 "algorithm/Math/Convolution/discrete_fourier_transform.hpp"



/**
 * @brief Discrete Fourier Transform（離散フーリエ変換）
 * @docs docs/Math/Convolution/discrete_fourier_transform.md
 */

#include <algorithm>
#line 11 "algorithm/Math/Convolution/discrete_fourier_transform.hpp"
#include <complex>
#include <type_traits>
#line 14 "algorithm/Math/Convolution/discrete_fourier_transform.hpp"

namespace algorithm {

namespace dft {

using D = double;

const D PI = std::acos(-1.0);

// Discrete Fourier Transform（離散フーリエ変換）. O(N^2).
void transform(std::vector<std::complex<D> > &a, bool inv = false) {
    if(a.size() == 0) return;
    const int n = a.size();
    std::vector<std::complex<D> > res(n, 0.0);
    D ang = 2 * PI / n;
    if(inv) ang = -ang;
    for(int i = 0; i < n; ++i) {
        D tmp = ang * i;
        for(int j = 0; j < n; ++j) res[i] += a[j] * std::polar<D>(1.0, tmp * j);
    }
    if(inv) {
        for(int i = 0; i < n; ++i) res[i] /= n;
    }
    a = res;
}

// 畳み込み．
// 数列a[], b[]に対して，c[i]=sum_{k=0}^{i} a[k]*b[i-k] となる数列c[]を求める．O(N^2).
template <typename Type>
std::vector<Type> convolve_naive(const std::vector<Type> &a, const std::vector<Type> &b) {
    const int n = a.size(), m = b.size();
    if(n == 0 or m == 0) return std::vector<Type>();
    std::vector<Type> res(n + m - 1, 0);
    for(int i = 0; i < n; ++i) {
        for(int j = 0; j < m; ++j) res[i + j] += a[i] * b[j];
    }
    return res;
}

// 畳み込み．
// 数列a[], b[]に対して，c[i]=sum_{k=0}^{i} a[k]*b[i-k] となる数列c[]を求める．O(N^2).
template <typename Type, typename std::enable_if_t<std::is_integral_v<Type>, bool> = false>
std::vector<Type> convolve(const std::vector<Type> &a, const std::vector<Type> &b) {
    if(a.size() == 0 or b.size() == 0) return std::vector<Type>();
    const int n = a.size() + b.size() - 1;
    std::vector<std::complex<D> > na(n, 0.0), nb(n, 0.0);
    std::copy(a.begin(), a.end(), na.begin());
    std::copy(b.begin(), b.end(), nb.begin());
    transform(na), transform(nb);
    for(int i = 0; i < n; ++i) na[i] *= nb[i];
    transform(na, true);
    std::vector<Type> res(n);
    for(int i = 0; i < n; ++i) res[i] = na[i].real() + (na[i].real() < 0.0 ? -0.5 : 0.5);
    return res;
}

// 畳み込み．
// 数列a[], b[]に対して，c[i]=sum_{k=0}^{i} a[k]*b[i-k] となる数列c[]を求める．O(N^2).
template <typename Type, typename std::enable_if_t<std::is_floating_point_v<Type>, bool> = false>
std::vector<Type> convolve(const std::vector<Type> &a, const std::vector<Type> &b) {
    if(a.size() == 0 or b.size() == 0) return std::vector<Type>();
    const int n = a.size() + b.size() - 1;
    std::vector<std::complex<D> > na(n, 0.0), nb(n, 0.0);
    std::copy(a.begin(), a.end(), na.begin());
    std::copy(b.begin(), b.end(), nb.begin());
    transform(na), transform(nb);
    for(int i = 0; i < n; ++i) na[i] *= nb[i];
    transform(na, true);
    std::vector<Type> res(n);
    for(int i = 0; i < n; ++i) res[i] = na[i].real();
    return res;
}

}  // namespace dft

}  // namespace algorithm


#line 10 "verify/aoj-ITP1_1_A-discrete_fourier_transform.test.cpp"

int main() {
    constexpr int t = 100;
    std::random_device seed;
    std::mt19937_64 rng(seed());

    constexpr double eps = 1e-6;
    std::uniform_int_distribution<int> uniform_n(1, 200);
    std::uniform_real_distribution<double> uniform(-1e3, 1e3);
    for(int i = 0; i < t; ++i) {
        const int n = uniform_n(rng), m = uniform_n(rng);
        std::vector<double> a(n), b(m);
        for(auto &elem : a) elem = uniform(rng);
        for(auto &elem : b) elem = uniform(rng);

        auto &&target = algorithm::dft::convolve(a, b);
        auto &&want = algorithm::dft::convolve_naive(a, b);

        assert(target.size() == size_t(n + m - 1));
        assert(want.size() == size_t(n + m - 1));
        for(int j = 0; j < n + m - 1; ++j) assert(std::abs(target[j] - want[j]) < eps);
    }

    std::cout << "Hello World" << std::endl;
}

Back to top page