// 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 <math/core/Complex.hpp>
#include <cmath>
#include <span>
#include <algorithm>
#include <type_traits>
#include <numbers>

#include <math/fft/fft_simd.hpp>

namespace sangi {

// Trait to detect whether a type is sangi::Complex (MSVC cannot handle
// partial specialization of concept-constrained classes, so detect via members)
namespace detail {
template<typename U>
constexpr bool is_sangi_complex_v = requires(U u) {
    u.re;
    u.im;
} && !std::is_same_v<U, std::complex<double>>
  && !std::is_same_v<U, std::complex<float>>
  && !std::is_same_v<U, std::complex<long double>>;

// For twiddle computation — branch with if constexpr between standard floating point and multi-precision (sangi::Float, etc.)
template<typename Real>
inline Real fft_pi() {
    if constexpr (std::is_floating_point_v<Real>) {
        return std::numbers::pi_v<Real>;
    } else {
        // Expects sangi::Float: Real::pi(prec) / Real::defaultPrecision()
        return Real::pi(Real::defaultPrecision());
    }
}
template<typename Real>
inline Real fft_cos(const Real& x) {
    if constexpr (std::is_floating_point_v<Real>) {
        return std::cos(x);
    } else {
        return cos(x, Real::defaultPrecision());
    }
}
template<typename Real>
inline Real fft_sin(const Real& x) {
    if constexpr (std::is_floating_point_v<Real>) {
        return std::sin(x);
    } else {
        return sin(x, Real::defaultPrecision());
    }
}
} // namespace detail

// FFT normalization mode
enum class FFTDivMode {
    None,       // no normalization
    Forward,    // normalize on the forward transform (1/N)
    Inverse,    // normalize on the inverse transform (1/N) — default
    Both        // normalize on both (1/sqrt(N))
};

/**
 * @brief Fast Fourier Transform (FFT) class
 *
 * This class provides a mixed-radix (radix 2, 3, 5) FFT.
 *
 * @tparam T value type (typically a complex type, or ModularInt)
 * @tparam Real real type (the Real type when T=Complex<Real>)
 */
template<typename T, typename Real = double>
class FFT {
private:
    // Type traits (complex type used for twiddle factors, etc.)
    using complex_type = Complex<Real>;

    // Detect whether the type is sangi::Complex
    template<typename U>
    static constexpr bool is_complex = detail::is_sangi_complex_v<U>;

    // Twiddle computation uses the template functions in namespace detail (branches to
    // std::cos/sin if Real is standard floating point such as float/double/long double,
    // or to Real::pi(prec) + cos/sin(x, prec) for sangi::Float, etc.)
    static Real fft_pi()                  { return detail::fft_pi<Real>(); }
    static Real fft_cos(const Real& x)    { return detail::fft_cos<Real>(x); }
    static Real fft_sin(const Real& x)    { return detail::fft_sin<Real>(x); }

    // Compute a point on the unit circle
    static complex_type omega(int k, int n) {
        const Real theta = Real(2) * fft_pi() * static_cast<Real>(k) / static_cast<Real>(n);
        return complex_type(fft_cos(theta), fft_sin(theta));
    }
    
    // Prime factorization
    static std::vector<int> factorize(int n) {
        std::vector<int> factors;

        // Extract factors of 5
        while (n % 5 == 0) {
            factors.push_back(5);
            n /= 5;
        }

        // Extract factors of 3
        while (n % 3 == 0) {
            factors.push_back(3);
            n /= 3;
        }

        // Extract factors of 2
        while (n % 2 == 0) {
            factors.push_back(2);
            n /= 2;
        }

        // Remaining prime factor (if any)
        if (n > 1) {
            factors.push_back(n);
        }

        // Sort factors in ascending order
        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);
        }
    }
    
    // Apply bit-reversal ordering
    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]);
            }
            
            // Bit-reversal index update
            int mask = n;
            do {
                mask >>= 1;
                j ^= mask;
            } while (mask && (j & mask) == 0);
        }
    }
    
    // Bluestein FFT: computes the DFT of arbitrary length N via radix-2 convolution (Chirp-Z transform)
    //
    // Using the identity nk = (n² + k² - (k-n)²) / 2,
    // the DFT is reduced to a circular convolution of length M ≥ 2N-1 (a power of 2).
    static void bluestein_fft(std::span<T> data, bool inverse = false) {
        const int N = static_cast<int>(data.size());
        if (N <= 1) return;

        // Smallest power of 2 with M ≥ 2N-1
        int M = 1;
        while (M < 2 * N - 1) M <<= 1;

        const Real pi = fft_pi();
        // Forward: exp(-πi·n²/N), inverse: exp(+πi·n²/N)
        const Real sign = inverse ? Real(1) : Real(-1);

        // Chirp sequence: 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(fft_cos(angle), fft_sin(angle));
        }

        // a[n] = x[n] · chirp[n], zero-padded to length M
        std::vector<complex_type> a(M, complex_type(Real(0), Real(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) arranged cyclically
        // b[0]   = conj(chirp[0])
        // b[n]   = conj(chirp[n])     (n = 1..N-1)
        // b[M-n] = conj(chirp[n])     (n = 1..N-1, wrap of negative indices)
        std::vector<complex_type> b(M, complex_type(Real(0), Real(0)));
        b[0] = conj(chirp[0]);
        for (int n = 1; n < N; ++n) {
            b[n] = conj(chirp[n]);
            b[M - n] = conj(chirp[n]);
        }

        // M is a power of 2 → FFT via the radix-2 path
        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());
            }
        }

        // Inverse: 1/N normalization
        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 Check whether a size is computable by FFT
     * @param n input size
     * @return true if computable
     */
    static bool is_valid_size(int n) {
        if (n <= 0) return false;

        // Check whether it decomposes into a product of powers of 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 Round up to a valid FFT size
     * @param n minimum required size
     * @return smallest valid FFT size
     */
    static int next_valid_size(int n) {
        while (!is_valid_size(n)) {
            ++n;
        }
        return n;
    }
    
    /**
     * @brief Forward transform (FFT)
     * @param data input/output data array
     */
    static void fft(std::span<T> data) {
        const int n = static_cast<int>(data.size());
        if (n <= 1) return;

        // Special optimization when the size is a power of 2
        if ((n & (n - 1)) == 0) {
            // For float/double, dispatch to the SIMD-optimized path
            if constexpr (std::is_same_v<T, Complex<float>> || std::is_same_v<T, Complex<double>>) {
                simd_fft::fft(data.data(), n);
            } else {
                // Non float/double (ModularInt, etc.): scalar radix-2
                apply_bit_reversal(data);

                for (int len = 2; len <= n; len <<= 1) {
                    const int half_len = len / 2;
                    const Real angle = Real(-2) * fft_pi() / static_cast<Real>(len);

                    for (int i = 0; i < n; i += len) {
                        complex_type w(Real(1), Real(0));
                        complex_type wn(fft_cos(angle), fft_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)) {
            // Mixed-radix FFT (2,3,5)
            std::vector<int> factors = factorize(n);
            mixed_radix_fft_impl(data, factors, false);
        } else {
            // Arbitrary length: Bluestein/Chirp-Z
            bluestein_fft(data, false);
        }
    }

    /**
     * @brief Inverse transform (IFFT)
     * @param data input/output data array
     */
    static void ifft(std::span<T> data) {
        const int n = static_cast<int>(data.size());
        if (n <= 1) return;

        // General approach: take the complex conjugate, run the FFT, then divide by N
        if constexpr (is_complex<T>) {
            // For complex types
            for (int i = 0; i < n; ++i) {
                data[i] = conj(data[i]);
            }
            
            fft(data);
            
            Real inv_n = Real(1) / static_cast<Real>(n);
            for (int i = 0; i < n; ++i) {
                data[i] = conj(data[i]) * inv_n;
            }
        } else {
            // For non-complex types (ModularInt, etc.)
            // Special optimization when the size is a power of 2
            if ((n & (n - 1)) == 0) {
                // Bit reversal is the same as in the ordinary FFT
                apply_bit_reversal(data);
                
                for (int len = 2; len <= n; len <<= 1) {
                    const int half_len = len / 2;
                    const Real angle = Real(2) * fft_pi() / static_cast<Real>(len); // inverse: e^{+i2πkn/N}
                    
                    for (int i = 0; i < n; i += len) {
                        complex_type w(Real(1), Real(0));
                        complex_type wn(fft_cos(angle), fft_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;
                        }
                    }
                }
                
                // Divide by 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)) {
                // Mixed-radix FFT (with the inverse flag)
                std::vector<int> factors = factorize(n);
                mixed_radix_fft_impl(data, factors, true);
            } else {
                // Arbitrary length: Bluestein/Chirp-Z (inverse)
                bluestein_fft(data, true);
            }
        }
    }
    
    /**
     * @brief Polynomial multiplication (convolution)
     * @param a input polynomial 1
     * @param b input polynomial 2
     * @return the product polynomial
     */
    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;

        // Compute a size suitable for FFT
        int fft_size = 1;
        while (fft_size < n) {
            fft_size *= 2;
        }

        // Extend the arrays
        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];

        // Run the FFT
        fft(std::span<T>(fa));
        fft(std::span<T>(fb));

        // Multiply in the spectral domain
        for (int i = 0; i < fft_size; ++i) {
            fa[i] = fa[i] * fb[i];
        }

        // Inverse FFT
        ifft(std::span<T>(fa));

        // Truncate the result
        fa.resize(n);
        return fa;
    }
    
    /**
     * @brief Number Theoretic Transform (NTT) - FFT for ModularInt
     * @param data input/output data
     * @param primitive_root 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");
        }

        // Reorder into bit-reversal order
        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 Inverse Number Theoretic Transform (INTT) - inverse FFT for ModularInt
     * @param data input/output data
     * @param primitive_root 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");
        }

        // Reorder into bit-reversal order
        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;
                }
            }
        }
        
        // Divide by N
        ModularInt<P> inv_n = ModularInt<P>(n).inverse();
        for (int i = 0; i < n; ++i) {
            data[i] = data[i] * inv_n;
        }
    }
    
    /**
     * @brief Polynomial multiplication via NTT (for ModularInt)
     * @param a input polynomial 1
     * @param b input polynomial 2
     * @param primitive_root primitive root
     * @return the product polynomial
     */
    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;
        
        // Compute a size suitable for NTT (a power of 2)
        int ntt_size = 1;
        while (ntt_size < n) {
            ntt_size *= 2;
        }

        // Extend the arrays
        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];

        // Run the NTT
        ntt(std::span<ModularInt<P>>(fa), primitive_root);
        ntt(std::span<ModularInt<P>>(fb), primitive_root);

        // Multiply in the spectral domain
        for (int i = 0; i < ntt_size; ++i) {
            fa[i] = fa[i] * fb[i];
        }

        // Inverse NTT
        intt(std::span<ModularInt<P>>(fa), primitive_root);

        // Truncate the result
        fa.resize(n);
        return fa;
    }
};

// Specialization for real types (real FFT)
template<typename Real>
class RealFFT {
private:
    using complex_type = Complex<Real>;

    // Twiddle computation helpers (same as in the FFT class; see detail::fft_pi/cos/sin for details)
    static Real fft_pi()                  { return detail::fft_pi<Real>(); }
    static Real fft_cos(const Real& x)    { return detail::fft_cos<Real>(x); }
    static Real fft_sin(const Real& x)    { return detail::fft_sin<Real>(x); }

public:
    /**
     * @brief Fast Fourier Transform for real data
     * @param data real data
     * @return complex spectrum (size is data.size()/2+1)
     */
    static std::vector<Complex<Real>> fft(std::span<const Real> data) {
        const int n = static_cast<int>(data.size());

        // Even N: Packing (N/2-point complex FFT + post-processing) — about 2x faster
        if (n >= 2 && (n & 1) == 0) {
            return fft_packing(data, n);
        }

        // Odd N: full-size complex FFT fallback
        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: pack N reals into N/2 complex numbers and 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-point complex FFT
        FFT<complex_type, Real>::fft(std::span<complex_type>(z));

        // Unpack: twiddle factors computed by recurrence (only 2 cos/sin calls)
        std::vector<complex_type> result(half + 1);

        // k=0, k=N/2 special cases
        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 recurrence: wk *= w_step
        const Real angle_step = Real(-2) * fft_pi() / static_cast<Real>(n);
        complex_type w_step(fft_cos(angle_step), fft_sin(angle_step));
        complex_type wk = w_step;

        for (int k = 1; k < half; ++k) {
            complex_type zk = z[k];
            complex_type zmk_conj = 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 Inverse transform from a complex spectrum to real data
     * @param spectrum complex spectrum (size is half the output size + 1)
     * @param output_size size of the output real data
     * @return real data
     */
    static std::vector<Real> ifft(std::span<const 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");

        // Even N: Packing inverse transform
        if (n >= 2 && (n & 1) == 0) {
            return ifft_packing(spectrum, n);
        }

        // Odd N: fallback
        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] = 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 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 recurrence: W^{-k} = exp(+2πik/N)
        const Real angle_step = Real(2) * fft_pi() / static_cast<Real>(n);
        complex_type w_step(fft_cos(angle_step), fft_sin(angle_step));
        complex_type wk_inv = w_step;

        for (int k = 1; k < half; ++k) {
            complex_type A = X[k];
            complex_type B = 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;
    }

    // Old ifft_packing (kept for reference, unused)
    [[maybe_unused]] static std::vector<Real> ifft_packing_old(std::span<const Complex<Real>> X, int n) {
        const int half = n / 2;
        const Real pi = fft_pi();

        // Decompose X[k] = Xe[k] + W^k · Xo[k] to recover 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)
        // However, since conjugate symmetry gives X[N-k] = conj(X[k]),
        // the conjugate of X[N-k] (k < N/2) = X[k] ... no,
        // X is the DFT of a real signal, so X[N-k] = conj(X[k]).
        // Therefore conj(X[N-k]) = X[k].
        // Xe[k] = (X[k] + X[k]) / 2 = X[k]? That is wrong.
        //
        // Correctly: the input X[k] is the spectrum of a real signal, and X[N-k] is recovered from the full length of X:
        // for k < N/2, X[N-k] = conj(X[k]) (conjugate symmetry)
        // hence conj(X[N-k]) = X[k]
        //
        // Preprocessing: recover Z[k] inversely
        // in the forward: X[k] = Xe[k] + W^k · Xo[k] where Z[k] = Xe[k] + i·Xo[k]
        // that is: Xe[k] = Re(Z[k]) + i·0, Xo[k] = Im(Z[k]) + i·0 ... no,
        // Xe, Xo are complex in general.
        //
        // Inverse computation: separate by multiplying by W^{-k}
        // X[k] = Xe[k] + W^k · Xo[k]
        // use conj(X[N/2-k]) (from the periodicity of 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])
        // ... this is complicated. Use the direct formula.

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

        // k=0: Z[0] = complex(X[0].real, X[N/2].real) / each component
        // 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(fft_cos(angle), fft_sin(angle));

            // Recover Xe, Xo from X[k] and X[N-k]=conj(X[k])
            complex_type xk = X[k];
            complex_type xnk = conj(X[k]);  // X[N-k] = conj(X[k]) for real signal

            complex_type xe = (xk + 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 = conj(X[mk]);

            complex_type xek = (A + B) / Real(2);
            const Real ang_k = Real(-2) * pi * static_cast<Real>(k) / static_cast<Real>(n);
            complex_type xok = (A - B) / (Real(2) * complex_type(fft_cos(ang_k), fft_sin(ang_k)));
            z[k] = xek + complex_type(0, 1) * xok;
        }

        // N/2-point 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 Convolution for real data
     * @param a input array 1
     * @param b input array 2
     * @return convolution result
     */
    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;
        
        // Compute a size suitable for FFT
        int fft_size = 1;
        while (fft_size < n) {
            fft_size *= 2;
        }

        // Extend the real arrays
        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];

        // Run the real FFT
        auto fa_spectrum = fft(std::span<const Real>(fa));
        auto fb_spectrum = fft(std::span<const Real>(fb));

        // Multiply in the spectral domain
        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];
        }

        // Inverse FFT
        std::vector<Real> result = ifft(std::span<const complex_type>(product_spectrum), fft_size);

        // Truncate the result
        result.resize(n);
        return result;
    }
};

/**
 * @brief Instance-based FFT wrapper
 *
 * Provides an API compatible with the old FFT.hpp.
 * Internally delegates to the static API of FFT<complex_type, Real> and RealFFT<Real>.
 *
 * @tparam Real real type (double, float)
 */
template<typename Real = double>
class FFTEngine {
    using complex_type = Complex<Real>;

    int m_nfft;
    FFTDivMode m_divMode;

    // Apply normalization
    // Note: the underlying FFT::ifft already includes 1/N normalization.
    // FFTEngine applies an additional correction on top of that.
    //
    // Forward: FFT::fft is unnormalized → scale as needed
    // Inverse: FFT::ifft is already 1/N normalized → correct according to DivMode
    //
    //            Forward correction   Inverse correction
    // None       1                    N (cancels the 1/N)
    // Forward    1/N                  N (cancels the 1/N)
    // Inverse    1                    1 (keeps the built-in 1/N)
    // Both       1/sqrt(N)            sqrt(N) (corrects 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 (keeps the built-in 1/N)
        }

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

public:
    /**
     * @brief Constructor
     * @param nfft FFT size (a power of 2)
     * @param mode normalization 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 Complex FFT forward transform */
    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 Complex FFT inverse transform */
    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 Real FFT forward transform (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 Real FFT inverse transform (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