// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// fft.hpp
#ifndef SANGI_FFT_HPP
#define SANGI_FFT_HPP

#include <math/core/common.hpp>
#include <math/core/modular/ModularInt.hpp>

#include <vector>
#include <complex>
#include <cmath>
#include <span>
#include <algorithm>
#include <type_traits>
#include <numbers>

#include <math/fft/fft_simd.hpp>

namespace sangi {

// FFT 正規化モード
enum class FFTDivMode {
    None,       // 正規化なし
    Forward,    // 順変換で正規化 (1/N)
    Inverse,    // 逆変換で正規化 (1/N) — デフォルト
    Both        // 両方で正規化 (1/sqrt(N))
};

/**
 * @brief 高速フーリエ変換 (FFT) のクラス
 * 
 * このクラスは混合基数（radix 2, 3, 5）FFTを提供します。
 * 
 * @tparam T 値の型（通常は複素数型、あるいはModularInt）
 * @tparam Real 実数型（T=std::complex<Real>の場合のReal型）
 */
template<typename T, typename Real = double>
class FFT {
private:
    // 型特性
    using complex_type = std::conditional_t<std::is_same_v<T, std::complex<Real>>, T, std::complex<Real>>;
    
    // 型が std::complex かどうかを判定
    template<typename U> struct is_complex_t : std::false_type {};
    template<typename R> struct is_complex_t<std::complex<R>> : std::true_type {};
    template<typename U> static constexpr bool is_complex = is_complex_t<U>::value;
    
    // 単位円上の点を計算
    static complex_type omega(int k, int n) {
        const Real theta = 2.0 * std::numbers::pi_v<Real> * static_cast<Real>(k) / static_cast<Real>(n);
        return complex_type(std::cos(theta), std::sin(theta));
    }
    
    // 素因数分解
    static std::vector<int> factorize(int n) {
        std::vector<int> factors;
        
        // 5の因数を抽出
        while (n % 5 == 0) {
            factors.push_back(5);
            n /= 5;
        }
        
        // 3の因数を抽出
        while (n % 3 == 0) {
            factors.push_back(3);
            n /= 3;
        }
        
        // 2の因数を抽出
        while (n % 2 == 0) {
            factors.push_back(2);
            n /= 2;
        }
        
        // 残りの素因数（ある場合）
        if (n > 1) {
            factors.push_back(n);
        }
        
        // 因数を小さい順にソート
        std::sort(factors.begin(), factors.end());
        
        return factors;
    }
    
    // Radix-2 butterfly with twiddle factors (Cooley-Tukey DIT)
    // stride encodes twiddle step; sign(stride) encodes forward(+) / inverse(-)
    static void butterfly_radix2(std::span<T> data, int stride, int offset, int m) {
        const int n = static_cast<int>(data.size());
        const int twiddle_stride = std::abs(stride);
        const int sign = (stride < 0) ? -1 : 1;

        for (int q = 0; q < m; ++q) {
            const int idx0 = offset + q;
            const int idx1 = idx0 + m;

            const long long exp1 = static_cast<long long>(sign) * q * twiddle_stride;
            complex_type tw1 = omega(static_cast<int>(exp1 % n + n) % n, n);

            complex_type x0 = static_cast<complex_type>(data[idx0]);
            complex_type x1 = static_cast<complex_type>(data[idx1]) * tw1;

            data[idx0] = static_cast<T>(x0 + x1);
            data[idx1] = static_cast<T>(x0 - x1);
        }
    }

    // Radix-3 butterfly with twiddle factors
    static void butterfly_radix3(std::span<T> data, int stride, int offset, int m) {
        const int n = static_cast<int>(data.size());
        const int twiddle_stride = std::abs(stride);
        const int sign = (stride < 0) ? -1 : 1;

        const complex_type w1 = omega(sign * 1, 3);
        const complex_type w2 = omega(sign * 2, 3);

        for (int q = 0; q < m; ++q) {
            const int idx0 = offset + q;
            const int idx1 = idx0 + m;
            const int idx2 = idx1 + m;

            auto tw = [&](int branch) -> complex_type {
                long long exp = static_cast<long long>(sign) * branch * q * twiddle_stride;
                return omega(static_cast<int>((exp % n + n) % n), n);
            };

            complex_type x0 = static_cast<complex_type>(data[idx0]);
            complex_type x1 = static_cast<complex_type>(data[idx1]) * tw(1);
            complex_type x2 = static_cast<complex_type>(data[idx2]) * tw(2);

            data[idx0] = static_cast<T>(x0 + x1 + x2);
            data[idx1] = static_cast<T>(x0 + w1 * x1 + w2 * x2);
            data[idx2] = static_cast<T>(x0 + w2 * x1 + w1 * x2);
        }
    }

    // Radix-5 butterfly with twiddle factors
    static void butterfly_radix5(std::span<T> data, int stride, int offset, int m) {
        const int n = static_cast<int>(data.size());
        const int twiddle_stride = std::abs(stride);
        const int sign = (stride < 0) ? -1 : 1;

        const complex_type w1 = omega(sign * 1, 5);
        const complex_type w2 = omega(sign * 2, 5);
        const complex_type w3 = omega(sign * 3, 5);
        const complex_type w4 = omega(sign * 4, 5);

        for (int q = 0; q < m; ++q) {
            const int idx0 = offset + q;
            const int idx1 = idx0 + m;
            const int idx2 = idx1 + m;
            const int idx3 = idx2 + m;
            const int idx4 = idx3 + m;

            auto tw = [&](int branch) -> complex_type {
                long long exp = static_cast<long long>(sign) * branch * q * twiddle_stride;
                return omega(static_cast<int>((exp % n + n) % n), n);
            };

            complex_type x0 = static_cast<complex_type>(data[idx0]);
            complex_type x1 = static_cast<complex_type>(data[idx1]) * tw(1);
            complex_type x2 = static_cast<complex_type>(data[idx2]) * tw(2);
            complex_type x3 = static_cast<complex_type>(data[idx3]) * tw(3);
            complex_type x4 = static_cast<complex_type>(data[idx4]) * tw(4);

            data[idx0] = static_cast<T>(x0 + x1 + x2 + x3 + x4);
            data[idx1] = static_cast<T>(x0 + w1 * x1 + w2 * x2 + w3 * x3 + w4 * x4);
            data[idx2] = static_cast<T>(x0 + w2 * x1 + w4 * x2 + w1 * x3 + w3 * x4);
            data[idx3] = static_cast<T>(x0 + w3 * x1 + w1 * x2 + w4 * x3 + w2 * x4);
            data[idx4] = static_cast<T>(x0 + w4 * x1 + w3 * x2 + w2 * x3 + w1 * x4);
        }
    }
    
    // ビット反転順序を適用
    template<typename U = T>
    static void apply_bit_reversal(std::span<U> data) {
        int n = static_cast<int>(data.size());
        for (int i = 0, j = 0; i < n; ++i) {
            if (i < j) {
                std::swap(data[i], data[j]);
            }
            
            // ビット反転のインデックス更新
            int mask = n;
            do {
                mask >>= 1;
                j ^= mask;
            } while (mask && (j & mask) == 0);
        }
    }
    
    // Bluestein FFT: 任意長 N の DFT を radix-2 畳み込みで計算 (Chirp-Z 変換)
    //
    // nk = (n² + k² - (k-n)²) / 2 の恒等式を利用して
    // DFT を長さ M ≥ 2N-1 (2のべき乗) の巡回畳み込みに帰着する。
    static void bluestein_fft(std::span<T> data, bool inverse = false) {
        const int N = static_cast<int>(data.size());
        if (N <= 1) return;

        // M ≥ 2N-1 の最小 2 のべき乗
        int M = 1;
        while (M < 2 * N - 1) M <<= 1;

        const Real pi = std::numbers::pi_v<Real>;
        // 順変換: exp(-πi·n²/N), 逆変換: exp(+πi·n²/N)
        const Real sign = inverse ? Real(1) : Real(-1);

        // チャープ列: chirp[n] = exp(sign·πi·n²/N)
        std::vector<complex_type> chirp(N);
        for (int n = 0; n < N; ++n) {
            Real angle = sign * pi
                * static_cast<Real>(static_cast<long long>(n) * n)
                / static_cast<Real>(N);
            chirp[n] = complex_type(std::cos(angle), std::sin(angle));
        }

        // a[n] = x[n] · chirp[n], ゼロパディングして長さ M
        std::vector<complex_type> a(M, complex_type(0, 0));
        for (int n = 0; n < N; ++n) {
            if constexpr (is_complex<T>) {
                a[n] = static_cast<complex_type>(data[n]) * chirp[n];
            } else {
                a[n] = complex_type(static_cast<Real>(data[n]), Real(0)) * chirp[n];
            }
        }

        // b: conj(chirp) を巡回配置
        // b[0]   = conj(chirp[0])
        // b[n]   = conj(chirp[n])     (n = 1..N-1)
        // b[M-n] = conj(chirp[n])     (n = 1..N-1, 負インデックスのラップ)
        std::vector<complex_type> b(M, complex_type(0, 0));
        b[0] = std::conj(chirp[0]);
        for (int n = 1; n < N; ++n) {
            b[n] = std::conj(chirp[n]);
            b[M - n] = std::conj(chirp[n]);
        }

        // M は 2 のべき乗 → radix-2 パスで FFT
        FFT<complex_type, Real>::fft(std::span<complex_type>(a));
        FFT<complex_type, Real>::fft(std::span<complex_type>(b));

        for (int i = 0; i < M; ++i) a[i] *= b[i];

        FFT<complex_type, Real>::ifft(std::span<complex_type>(a));

        // X[k] = chirp[k] · (a ★ b)[k]
        for (int k = 0; k < N; ++k) {
            complex_type val = chirp[k] * a[k];
            if constexpr (is_complex<T>) {
                data[k] = val;
            } else {
                data[k] = static_cast<T>(val.real());
            }
        }

        // 逆変換: 1/N 正規化
        if (inverse) {
            Real inv_n = Real(1) / static_cast<Real>(N);
            for (int k = 0; k < N; ++k) {
                if constexpr (is_complex<T>) {
                    data[k] *= inv_n;
                } else {
                    data[k] = static_cast<T>(static_cast<Real>(data[k]) * inv_n);
                }
            }
        }
    }

    // Mixed-radix FFT (Cooley-Tukey DIT with digit-reversal permutation)
    // Factors: sorted list of radices (2, 3, 5). Product must equal data.size().
    static void mixed_radix_fft_impl(std::span<T> data, const std::vector<int>& factors, bool inverse = false) {
        const int n = static_cast<int>(data.size());
        if (n <= 1) return;

        int product = 1;
        for (int factor : factors) {
            if (factor != 2 && factor != 3 && factor != 5)
                throw MathError("Unsupported radix factor");
            product *= factor;
        }
        if (product != n)
            throw MathError("Input size must be a product of supported radices (2, 3, 5)");

        // Mixed-radix digit-reversal permutation
        {
            std::vector<T> permuted(n);
            std::vector<int> digits(factors.size());
            for (int i = 0; i < n; ++i) {
                int x = i;
                for (std::size_t d = 0; d < factors.size(); ++d) {
                    digits[d] = x % factors[d];
                    x /= factors[d];
                }
                int j = 0, place = 1;
                for (int d = static_cast<int>(factors.size()) - 1; d >= 0; --d) {
                    j += digits[d] * place;
                    place *= factors[d];
                }
                permuted[j] = data[i];
            }
            std::copy(permuted.begin(), permuted.end(), data.begin());
        }

        // Butterfly stages (iterate factors in reverse for DIT)
        int m = 1;
        for (auto it = factors.rbegin(); it != factors.rend(); ++it) {
            const int factor = *it;
            const int block_size = m * factor;
            const int stride = n / block_size;
            // SIMD radix-2 uses e^{-2πi/N} for forward, so mixed-radix must match
            const int signed_stride = inverse ? stride : -stride;

            for (int offset = 0; offset < n; offset += block_size) {
                if (factor == 2) {
                    butterfly_radix2(data, signed_stride, offset, m);
                } else if (factor == 3) {
                    butterfly_radix3(data, signed_stride, offset, m);
                } else {
                    butterfly_radix5(data, signed_stride, offset, m);
                }
            }
            m = block_size;
        }

        // Inverse: scale by 1/N
        if (inverse) {
            const Real inv_n = Real(1) / static_cast<Real>(n);
            for (int i = 0; i < n; ++i) {
                if constexpr (is_complex<T>) {
                    data[i] *= inv_n;
                } else {
                    data[i] = static_cast<T>(static_cast<Real>(data[i]) * inv_n);
                }
            }
        }
    }
    
public:
    /**
     * @brief FFT計算が可能なサイズかを確認
     * @param n 入力サイズ
     * @return 計算可能な場合はtrue
     */
    static bool is_valid_size(int n) {
        if (n <= 0) return false;
        
        // 2, 3, 5の累乗の積に分解できるか確認
        while (n % 2 == 0) n /= 2;
        while (n % 3 == 0) n /= 3;
        while (n % 5 == 0) n /= 5;
        
        return n == 1;
    }
    
    /**
     * @brief 有効なFFTサイズに切り上げ
     * @param n 必要最小サイズ
     * @return 最小の有効なFFTサイズ
     */
    static int next_valid_size(int n) {
        while (!is_valid_size(n)) {
            ++n;
        }
        return n;
    }
    
    /**
     * @brief 順変換 (FFT)
     * @param data 入力/出力データ配列
     */
    static void fft(std::span<T> data) {
        const int n = static_cast<int>(data.size());
        if (n <= 1) return;
        
        // サイズが2の累乗の場合は特殊最適化
        if ((n & (n - 1)) == 0) {
            // float/double の場合は SIMD 最適化パスにディスパッチ
            if constexpr (std::is_same_v<T, std::complex<float>> || std::is_same_v<T, std::complex<double>>) {
                simd_fft::fft(data.data(), n);
            } else {
                // 非 float/double (ModularInt 等): スカラー radix-2
                apply_bit_reversal(data);

                for (int len = 2; len <= n; len <<= 1) {
                    const int half_len = len / 2;
                    const Real angle = 2.0 * std::numbers::pi_v<Real> / static_cast<Real>(len);

                    for (int i = 0; i < n; i += len) {
                        complex_type w(1.0, 0.0);
                        complex_type wn(std::cos(angle), std::sin(angle));

                        for (int j = 0; j < half_len; ++j) {
                            T u = data[i + j];
                            T v = data[i + j + half_len] * w;

                            data[i + j] = u + v;
                            data[i + j + half_len] = u - v;

                            w *= wn;
                        }
                    }
                }
            }
        } else if (is_valid_size(n)) {
            // 混合基数FFT (2,3,5)
            std::vector<int> factors = factorize(n);
            mixed_radix_fft_impl(data, factors, false);
        } else {
            // 任意長: Bluestein/Chirp-Z
            bluestein_fft(data, false);
        }
    }
    
    /**
     * @brief 逆変換 (IFFT)
     * @param data 入力/出力データ配列
     */
    static void ifft(std::span<T> data) {
        const int n = static_cast<int>(data.size());
        if (n <= 1) return;
        
        // 一般的なアプローチ: 複素共役を取り、FFTを実行し、Nで割る
        if constexpr (is_complex<T>) {
            // 複素数型の場合
            for (int i = 0; i < n; ++i) {
                data[i] = std::conj(data[i]);
            }
            
            fft(data);
            
            Real inv_n = 1.0 / static_cast<Real>(n);
            for (int i = 0; i < n; ++i) {
                data[i] = std::conj(data[i]) * inv_n;
            }
        } else {
            // 非複素数型の場合（ModularIntなど）
            // サイズが2の累乗の場合は特殊最適化
            if ((n & (n - 1)) == 0) {
                // ビット反転は通常のFFTと同じ
                apply_bit_reversal(data);
                
                for (int len = 2; len <= n; len <<= 1) {
                    const int half_len = len / 2;
                    const Real angle = -2.0 * std::numbers::pi_v<Real> / static_cast<Real>(len); // 負の角度
                    
                    for (int i = 0; i < n; i += len) {
                        complex_type w(1.0, 0.0);
                        complex_type wn(std::cos(angle), std::sin(angle));
                        
                        for (int j = 0; j < half_len; ++j) {
                            T u = data[i + j];
                            T v = data[i + j + half_len] * w;
                            
                            data[i + j] = u + v;
                            data[i + j + half_len] = u - v;
                            
                            w *= wn;
                        }
                    }
                }
                
                // Nで割る
                T inv_n = static_cast<T>(1) / static_cast<T>(n);
                for (int i = 0; i < n; ++i) {
                    data[i] = data[i] * inv_n;
                }
            } else if (is_valid_size(n)) {
                // 混合基数FFT（逆変換フラグ付き）
                std::vector<int> factors = factorize(n);
                mixed_radix_fft_impl(data, factors, true);
            } else {
                // 任意長: Bluestein/Chirp-Z (逆変換)
                bluestein_fft(data, true);
            }
        }
    }
    
    /**
     * @brief 多項式乗算（畳み込み）
     * @param a 入力多項式1
     * @param b 入力多項式2
     * @return 乗算結果の多項式
     */
    static std::vector<T> convolve(std::span<const T> a, std::span<const T> b) {
        const int na = static_cast<int>(a.size());
        const int nb = static_cast<int>(b.size());
        const int n = na + nb - 1;
        
        // FFTに適したサイズを計算
        int fft_size = 1;
        while (fft_size < n) {
            fft_size *= 2;
        }
        
        // 配列を拡張
        std::vector<T> fa(fft_size, static_cast<T>(0));
        std::vector<T> fb(fft_size, static_cast<T>(0));
        
        for (int i = 0; i < na; ++i) fa[i] = a[i];
        for (int i = 0; i < nb; ++i) fb[i] = b[i];
        
        // FFTを実行
        fft(std::span<T>(fa));
        fft(std::span<T>(fb));
        
        // スペクトル領域で乗算
        for (int i = 0; i < fft_size; ++i) {
            fa[i] = fa[i] * fb[i];
        }
        
        // 逆FFT
        ifft(std::span<T>(fa));
        
        // 結果を切り詰める
        fa.resize(n);
        return fa;
    }
    
    /**
     * @brief 高速数論変換 (NTT) - ModularInt用のFFT
     * @param data 入力/出力データ
     * @param primitive_root 原始根
     */
    template<int P>
    static void ntt(std::span<ModularInt<P>> data, ModularInt<P> primitive_root) {
        const int n = static_cast<int>(data.size());
        if (n <= 1) return;
        
        if ((n & (n - 1)) != 0) {
            throw MathError("NTT only supports power of 2 sizes");
        }
        
        // ビット反転順序に並べ替え
        apply_bit_reversal(data);
        
        for (int len = 2; len <= n; len <<= 1) {
            const int half_len = len / 2;
            const ModularInt<P> wn = primitive_root.pow((P - 1) / len);
            
            for (int i = 0; i < n; i += len) {
                ModularInt<P> w(1);
                
                for (int j = 0; j < half_len; ++j) {
                    ModularInt<P> u = data[i + j];
                    ModularInt<P> v = data[i + j + half_len] * w;
                    
                    data[i + j] = u + v;
                    data[i + j + half_len] = u - v;
                    
                    w = w * wn;
                }
            }
        }
    }
    
    /**
     * @brief 逆高速数論変換 (INTT) - ModularInt用の逆FFT
     * @param data 入力/出力データ
     * @param primitive_root 原始根
     */
    template<int P>
    static void intt(std::span<ModularInt<P>> data, ModularInt<P> primitive_root) {
        const int n = static_cast<int>(data.size());
        if (n <= 1) return;
        
        if ((n & (n - 1)) != 0) {
            throw MathError("INTT only supports power of 2 sizes");
        }
        
        // ビット反転順序に並べ替え
        apply_bit_reversal(data);
        
        for (int len = 2; len <= n; len <<= 1) {
            const int half_len = len / 2;
            const ModularInt<P> wn = primitive_root.pow((P - 1) / len).inverse();
            
            for (int i = 0; i < n; i += len) {
                ModularInt<P> w(1);
                
                for (int j = 0; j < half_len; ++j) {
                    ModularInt<P> u = data[i + j];
                    ModularInt<P> v = data[i + j + half_len] * w;
                    
                    data[i + j] = u + v;
                    data[i + j + half_len] = u - v;
                    
                    w = w * wn;
                }
            }
        }
        
        // Nで割る
        ModularInt<P> inv_n = ModularInt<P>(n).inverse();
        for (int i = 0; i < n; ++i) {
            data[i] = data[i] * inv_n;
        }
    }
    
    /**
     * @brief 多項式乗算をNTTで実行（ModularInt用）
     * @param a 入力多項式1
     * @param b 入力多項式2
     * @param primitive_root 原始根
     * @return 乗算結果の多項式
     */
    template<int P>
    static std::vector<ModularInt<P>> convolve_ntt(
        std::span<const ModularInt<P>> a, 
        std::span<const ModularInt<P>> b,
        ModularInt<P> primitive_root) {
        
        const int na = static_cast<int>(a.size());
        const int nb = static_cast<int>(b.size());
        const int n = na + nb - 1;
        
        // NTTに適したサイズを計算（2の累乗）
        int ntt_size = 1;
        while (ntt_size < n) {
            ntt_size *= 2;
        }
        
        // 配列を拡張
        std::vector<ModularInt<P>> fa(ntt_size, ModularInt<P>(0));
        std::vector<ModularInt<P>> fb(ntt_size, ModularInt<P>(0));
        
        for (int i = 0; i < na; ++i) fa[i] = a[i];
        for (int i = 0; i < nb; ++i) fb[i] = b[i];
        
        // NTTを実行
        ntt(std::span<ModularInt<P>>(fa), primitive_root);
        ntt(std::span<ModularInt<P>>(fb), primitive_root);
        
        // スペクトル領域で乗算
        for (int i = 0; i < ntt_size; ++i) {
            fa[i] = fa[i] * fb[i];
        }
        
        // 逆NTT
        intt(std::span<ModularInt<P>>(fa), primitive_root);
        
        // 結果を切り詰める
        fa.resize(n);
        return fa;
    }
};

// 実数型向けの特殊化（実FFT）
template<typename Real>
class RealFFT {
private:
    using complex_type = std::complex<Real>;
    
public:
    /**
     * @brief 実数データに対する高速フーリエ変換
     * @param data 実数データ
     * @return 複素数スペクトル（サイズはdata.size()/2+1）
     */
    static std::vector<std::complex<Real>> fft(std::span<const Real> data) {
        const int n = static_cast<int>(data.size());

        // 偶数 N: Packing (N/2 点複素 FFT + 後処理) — 約 2x 高速
        if (n >= 2 && (n & 1) == 0) {
            return fft_packing(data, n);
        }

        // 奇数 N: フルサイズ複素 FFT フォールバック
        std::vector<complex_type> complex_data(n);
        for (int i = 0; i < n; ++i)
            complex_data[i] = complex_type(data[i], 0);
        FFT<complex_type, Real>::fft(std::span<complex_type>(complex_data));
        std::vector<complex_type> result(n / 2 + 1);
        for (int i = 0; i <= n / 2; ++i)
            result[i] = complex_data[i];
        return result;
    }

private:
    // Packing: N 個の実数を N/2 個の複素数に詰めて FFT
    // z[k] = x[2k] + i·x[2k+1], FFT(z, N/2), then unpack
    static std::vector<complex_type> fft_packing(std::span<const Real> data, int n) {
        const int half = n / 2;

        // Pack: z[k] = x[2k] + i·x[2k+1]
        std::vector<complex_type> z(half);
        for (int k = 0; k < half; ++k)
            z[k] = complex_type(data[2*k], data[2*k + 1]);

        // N/2 点複素 FFT
        FFT<complex_type, Real>::fft(std::span<complex_type>(z));

        // Unpack: twiddle 因子は漸化式で計算 (cos/sin 呼び出し 2 回のみ)
        std::vector<complex_type> result(half + 1);

        // k=0, k=N/2 特殊ケース
        result[0] = complex_type(z[0].real() + z[0].imag(), 0);
        result[half] = complex_type(z[0].real() - z[0].imag(), 0);

        // W_N = exp(-2πi/N), twiddle 漸化式: wk *= w_step
        const Real angle_step = Real(-2) * std::numbers::pi_v<Real> / n;
        complex_type w_step(std::cos(angle_step), std::sin(angle_step));
        complex_type wk = w_step;

        for (int k = 1; k < half; ++k) {
            complex_type zk = z[k];
            complex_type zmk_conj = std::conj(z[half - k]);
            complex_type xe = (zk + zmk_conj) * Real(0.5);
            complex_type diff = zk - zmk_conj;
            complex_type xo(diff.imag() * Real(0.5), -diff.real() * Real(0.5));
            result[k] = xe + wk * xo;
            wk *= w_step;
        }

        return result;
    }

public:
    
    /**
     * @brief 複素スペクトルから実数データへの逆変換
     * @param spectrum 複素スペクトル（サイズは出力サイズの半分+1）
     * @param output_size 出力実数データのサイズ
     * @return 実数データ
     */
    static std::vector<Real> ifft(std::span<const std::complex<Real>> spectrum, int output_size) {
        const int half_n_plus_1 = static_cast<int>(spectrum.size());
        const int n = output_size;

        if (half_n_plus_1 != n / 2 + 1)
            throw MathError("Spectrum size must be output_size/2+1");

        // 偶数 N: Packing 逆変換
        if (n >= 2 && (n & 1) == 0) {
            return ifft_packing(spectrum, n);
        }

        // 奇数 N: フォールバック
        std::vector<complex_type> complex_data(n);
        for (int i = 0; i < half_n_plus_1; ++i)
            complex_data[i] = spectrum[i];
        for (int i = 1; i < n - half_n_plus_1 + 1; ++i)
            complex_data[n - i] = std::conj(spectrum[i]);
        FFT<complex_type, Real>::ifft(std::span<complex_type>(complex_data));
        std::vector<Real> result(n);
        for (int i = 0; i < n; ++i)
            result[i] = complex_data[i].real();
        return result;
    }

private:
    // Inverse Packing: X[0..N/2] → Z[0..N/2-1] → IFFT → unpack
    // Inverse Packing: X[0..N/2] → Z[0..N/2-1] → IFFT → unpack
    // A=X[k], B=conj(X[N/2-k]) → Xe=(A+B)/2, W^k·Xo=(A-B)/2 → Z[k]=Xe+i·Xo
    static std::vector<Real> ifft_packing(std::span<const std::complex<Real>> X, int n) {
        const int half = n / 2;
        std::vector<complex_type> z(half);

        // k=0
        z[0] = complex_type(
            (X[0].real() + X[half].real()) / 2,
            (X[0].real() - X[half].real()) / 2
        );

        // twiddle 漸化式: W^{-k} = exp(+2πik/N)
        const Real angle_step = Real(2) * std::numbers::pi_v<Real> / n;
        complex_type w_step(std::cos(angle_step), std::sin(angle_step));
        complex_type wk_inv = w_step;

        for (int k = 1; k < half; ++k) {
            complex_type A = X[k];
            complex_type B = std::conj(X[half - k]);
            complex_type xe = (A + B) * Real(0.5);
            complex_type wk_xo = (A - B) * Real(0.5);  // W^k · Xo[k]
            complex_type xo = wk_inv * wk_xo;           // W^{-k} · W^k · Xo = Xo
            z[k] = xe + complex_type(-xo.imag(), xo.real()); // xe + i·xo
            wk_inv *= w_step;
        }

        FFT<complex_type, Real>::ifft(std::span<complex_type>(z));

        std::vector<Real> result(n);
        for (int k = 0; k < half; ++k) {
            result[2*k]     = z[k].real();
            result[2*k + 1] = z[k].imag();
        }
        return result;
    }

    // 旧 ifft_packing (参考用、未使用)
    [[maybe_unused]] static std::vector<Real> ifft_packing_old(std::span<const std::complex<Real>> X, int n) {
        const int half = n / 2;
        const Real pi = std::numbers::pi_v<Real>;

        // X[k] = Xe[k] + W^k · Xo[k] を分解して Z[k] = Xe[k] + i·Xo[k] を復元
        // Xe[k] = (X[k] + conj(X[N-k])) / 2
        // Xo[k] = (X[k] - conj(X[N-k])) / (2 · W^k)
        // ただし X[N-k] = conj(X[k]) が共役対称で成り立つので
        // X[N-k] (k < N/2) の共役 = X[k] ... いいえ、
        // X は実信号の DFT なので X[N-k] = conj(X[k])。
        // 従って conj(X[N-k]) = X[k]。
        // Xe[k] = (X[k] + X[k]) / 2 = X[k]? それは間違い。
        //
        // 正しくは: 入力の X[k] は実信号のスペクトルで、X[N-k] を X の全長から復元する:
        // k < N/2 のとき X[N-k] = conj(X[k]) (共役対称)
        // よって conj(X[N-k]) = X[k]
        //
        // 前処理: Z[k] を逆に求める
        // forward では: X[k] = Xe[k] + W^k · Xo[k] where Z[k] = Xe[k] + i·Xo[k]
        // つまり: Xe[k] = Re(Z[k]) + i·0, Xo[k] = Im(Z[k]) + i·0 ... いいえ、
        // Xe, Xo は一般に複素数。
        //
        // 逆算: W^{-k} を掛けて分離
        // X[k] = Xe[k] + W^k · Xo[k]
        // conj(X[N/2-k]) を使う (Z の周期性から):
        //   X_rev[k] := conj(X[N-k]) = Xe_conj_rev[k] + conj(W^{N-k}) · Xo_conj_rev[k]
        //   = conj(Xe[N/2-k]) + W^{k-N} · conj(Xo[N/2-k])
        // ... これは複雑。直接公式を使う。

        std::vector<complex_type> z(half);

        // k=0: Z[0] = complex(X[0].real, X[N/2].real) / 各成分
        // forward: X[0] = Re(Z[0]) + Im(Z[0]), X[N/2] = Re(Z[0]) - Im(Z[0])
        // inverse: Re(Z[0]) = (X[0] + X[N/2]) / 2, Im(Z[0]) = (X[0] - X[N/2]) / 2
        z[0] = complex_type(
            (X[0].real() + X[half].real()) / 2,
            (X[0].real() - X[half].real()) / 2
        );

        for (int k = 1; k < half; ++k) {
            // W_N^k = exp(-2πik/N) → W^{-k} = exp(+2πik/N)
            Real angle = 2 * pi * k / n;
            complex_type wk_inv(std::cos(angle), std::sin(angle));

            // X[k] と X[N-k]=conj(X[k]) から Xe, Xo を復元
            complex_type xk = X[k];
            complex_type xnk = std::conj(X[k]);  // X[N-k] = conj(X[k]) for real signal

            complex_type xe = (xk + std::conj(xnk)) / Real(2);  // = (X[k] + conj(conj(X[k]))) / 2 = X[k]
            // ... wait this simplifies to xe = xk. That's wrong.

            // Actually: X[N-k] for general k should use the full spectrum.
            // For 0 < k < N/2: X[N-k] is at index N-k > N/2, which is conj(X[k]).
            // But in our unpack formula, the "conj(X[N/2-k])" uses Z's period N/2.
            // Let me re-derive.

            // From forward: Z is the N/2-point FFT of z[n] = x[2n] + i·x[2n+1]
            // We have Z[k] and Z[N/2-k].
            // Xe[k] = (Z[k] + conj(Z[N/2-k])) / 2
            // diff = Z[k] - conj(Z[N/2-k])
            // Xo[k] = complex(diff.imag, -diff.real) / 2
            // X[k] = Xe[k] + W^k · Xo[k]
            //
            // To invert: we need Z[k] from X[k] and X[N-k]=conj(X[k]).
            // X[k] = Xe[k] + W^k · Xo[k]
            // X[N-k] = Xe[N/2-k'] ... this is complex. Use the pair approach:
            // X[k]         = Xe[k] + W^k  · Xo[k]
            // conj(X[N-k]) = Xe[k] - W^k  · Xo[k]   (from conjugate symmetry of Xe, Xo)
            // Wait, is Xe conjugate-symmetric? Xe[k] = DFT of x_even[n], which is real.
            // So Xe[N/2-k] = conj(Xe[k]). Similarly Xo[N/2-k] = conj(Xo[k]).

            // Pair equations:
            // X[k]             = Xe[k] + W^k · Xo[k]
            // conj(X[N-k])     = Xe[k] + conj(W^{N-k}) · Xo[k]
            //                   = Xe[k] + conj(W^{-k}) · Xo[k]
            //                   = Xe[k] + W^k · Xo[k]  (since |W|=1, conj(W^{-k}) = W^k)
            // ... that gives X[k] = conj(X[N-k]) which is just conjugate symmetry. Not useful.

            // I need to use Z[k] and Z[N/2-k] relationship, not X.
            // The correct inverse packing formula:
            // Given X[k] for k=0..N/2, recover Z[k] for k=0..N/2-1:
            //
            // A[k] = X[k]
            // B[k] = X[N/2 - k]  (or for k=0, B[0] = X[N/2])... hmm still complex.

            // Let me just use the direct formula. The forward computes X from Z.
            // For the inverse, I compute Z from X by inverting those formulas.
            // Since we stored result[k] = xe + wk*xo, and result[half-k] = ...,
            // I can set up a 2x2 system.

            // Actually, the simplest approach: compute Z[k] directly.
            // Z[k] = Xe[k] + i·Xo[k]
            // where Xe and Xo are real-signal DFTs.
            // From the forward post-processing:
            //   result[k] = Xe[k] + W^k · Xo[k]
            //   result[N-k] = conj(Xe[k] - W^k · Xo[k])  (using conj symmetry)
            // But N-k > N/2 so result[N-k] = conj(result[k'])...
            // I can also use result[N/2-k]:
            //   result[N/2-k] = Xe[N/2-k] + W^{N/2-k} · Xo[N/2-k]
            // Since Xe is periodic with period N/2... no wait, Xe is the DFT of even samples
            // with length N/2, so Xe has period N/2. So Xe[N/2-k] = Xe[-k] = conj(Xe[k]).
            // Similarly Xo[N/2-k] = conj(Xo[k]).
            // W^{N/2-k} = W^{N/2} · W^{-k} = -1 · conj(W^k) = -conj(W^k)
            //
            // So: result[N/2-k] = conj(Xe[k]) - conj(W^k) · conj(Xo[k])
            //                   = conj(Xe[k] - W^k · Xo[k])
            //                   = conj(Xe[k] - W^k · Xo[k])
            //
            // Let A = result[k] = Xe[k] + W^k · Xo[k]
            // Let B = conj(result[N/2-k]) = Xe[k] - W^k · Xo[k]
            //
            // Solving: Xe[k] = (A + B) / 2
            //          Xo[k] = (A - B) / (2 · W^k)
            //          Z[k] = Xe[k] + i · Xo[k]

            int mk = half - k;
            complex_type A = X[k];
            complex_type B = std::conj(X[mk]);

            complex_type xek = (A + B) / Real(2);
            complex_type xok = (A - B) / (Real(2) * complex_type(std::cos(-2*pi*k/n), std::sin(-2*pi*k/n)));
            z[k] = xek + complex_type(0, 1) * xok;
        }

        // N/2 点 IFFT
        FFT<complex_type, Real>::ifft(std::span<complex_type>(z));

        // Unpack: x[2k] = Re(z[k]), x[2k+1] = Im(z[k])
        std::vector<Real> result(n);
        for (int k = 0; k < half; ++k) {
            result[2*k]     = z[k].real();
            result[2*k + 1] = z[k].imag();
        }
        return result;
    }

public:
    
    /**
     * @brief 実数データに対する畳み込み
     * @param a 入力配列1
     * @param b 入力配列2
     * @return 畳み込み結果
     */
    static std::vector<Real> convolve(std::span<const Real> a, std::span<const Real> b) {
        const int na = static_cast<int>(a.size());
        const int nb = static_cast<int>(b.size());
        const int n = na + nb - 1;
        
        // FFTに適したサイズを計算
        int fft_size = 1;
        while (fft_size < n) {
            fft_size *= 2;
        }
        
        // 実数配列を拡張
        std::vector<Real> fa(fft_size, 0);
        std::vector<Real> fb(fft_size, 0);
        
        for (int i = 0; i < na; ++i) fa[i] = a[i];
        for (int i = 0; i < nb; ++i) fb[i] = b[i];
        
        // 実FFTを実行
        auto fa_spectrum = fft(std::span<const Real>(fa));
        auto fb_spectrum = fft(std::span<const Real>(fb));
        
        // スペクトル領域で乗算
        std::vector<complex_type> product_spectrum(fa_spectrum.size());
        for (size_t i = 0; i < fa_spectrum.size(); ++i) {
            product_spectrum[i] = fa_spectrum[i] * fb_spectrum[i];
        }
        
        // 逆FFT
        std::vector<Real> result = ifft(std::span<const complex_type>(product_spectrum), fft_size);
        
        // 結果を切り詰める
        result.resize(n);
        return result;
    }
};

/**
 * @brief インスタンスベースの FFT ラッパー
 *
 * 旧 FFT.hpp 互換の API を提供する。
 * 内部で FFT<complex_type, Real> および RealFFT<Real> の static API に委譲する。
 *
 * @tparam Real 実数型 (double, float)
 */
template<typename Real = double>
class FFTEngine {
    using complex_type = std::complex<Real>;

    int m_nfft;
    FFTDivMode m_divMode;

    // 正規化を適用
    // 注: 下層の FFT::ifft は既に 1/N 正規化を含む。
    // FFTEngine はその上で追加の補正を行う。
    //
    // Forward (順変換): FFT::fft は非正規化 → 必要に応じてスケーリング
    // Inverse (逆変換): FFT::ifft は 1/N 正規化済み → DivMode に応じて補正
    //
    //            Forward補正     Inverse補正
    // None       1               N (1/N を打ち消す)
    // Forward    1/N             N (1/N を打ち消す)
    // Inverse    1               1 (組み込み 1/N をそのまま)
    // Both       1/sqrt(N)       sqrt(N) (1/N → 1/sqrt(N) に補正)
    template<typename Container>
    void applyNormalization(Container& data, bool isForward) const {
        Real scale = Real(1);

        if (isForward) {
            if (m_divMode == FFTDivMode::Forward) {
                scale = Real(1) / static_cast<Real>(m_nfft);
            } else if (m_divMode == FFTDivMode::Both) {
                scale = Real(1) / std::sqrt(static_cast<Real>(m_nfft));
            }
        } else {
            if (m_divMode == FFTDivMode::None || m_divMode == FFTDivMode::Forward) {
                scale = static_cast<Real>(m_nfft);
            } else if (m_divMode == FFTDivMode::Both) {
                scale = std::sqrt(static_cast<Real>(m_nfft));
            }
            // Inverse: scale = 1 (組み込み 1/N がそのまま)
        }

        if (scale != Real(1)) {
            for (auto& x : data) {
                x *= scale;
            }
        }
    }

public:
    /**
     * @brief コンストラクタ
     * @param nfft FFT サイズ (2のべき乗)
     * @param mode 正規化モード
     */
    explicit FFTEngine(int nfft, FFTDivMode mode = FFTDivMode::Forward)
        : m_nfft(nfft), m_divMode(mode)
    {
        if (nfft <= 0 || (nfft & (nfft - 1)) != 0) {
            throw std::invalid_argument("FFT size must be a power of 2");
        }
    }

    int size() const { return m_nfft; }

    /** @brief 複素 FFT 順変換 */
    std::vector<complex_type> transform(const std::vector<complex_type>& input) {
        if (static_cast<int>(input.size()) != m_nfft) {
            throw std::invalid_argument("Input size must match FFT size");
        }
        std::vector<complex_type> output = input;
        FFT<complex_type, Real>::fft(std::span<complex_type>(output));
        applyNormalization(output, true);
        return output;
    }

    /** @brief 複素 FFT 逆変換 */
    std::vector<complex_type> inverse(const std::vector<complex_type>& input) {
        if (static_cast<int>(input.size()) != m_nfft) {
            throw std::invalid_argument("Input size must match FFT size");
        }
        std::vector<complex_type> output = input;
        FFT<complex_type, Real>::ifft(std::span<complex_type>(output));
        applyNormalization(output, false);
        return output;
    }

    /** @brief 実数 FFT 順変換 (N → N/2+1) */
    std::vector<complex_type> real_transform(const std::vector<Real>& input) {
        if (static_cast<int>(input.size()) != m_nfft) {
            throw std::invalid_argument("Input size must match FFT size");
        }
        auto output = RealFFT<Real>::fft(std::span<const Real>(input));
        applyNormalization(output, true);
        return output;
    }

    /** @brief 実数 FFT 逆変換 (N/2+1 → N) */
    std::vector<Real> real_inverse(const std::vector<complex_type>& input) {
        if (static_cast<int>(input.size()) != m_nfft / 2 + 1) {
            throw std::invalid_argument("Input size must be N/2+1 for real IFFT");
        }
        auto output = RealFFT<Real>::ifft(std::span<const complex_type>(input), m_nfft);
        applyNormalization(output, false);
        return output;
    }
};

} // namespace sangi

#endif // SANGI_FFT_HPP