/* Copyright (C) 2025 Joel Yliluoma — Material for https://youtu.be/C4tHxouarGc - https://iki.fi/bisqwit/ */
/* MIT license follows:
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
documentation files (the "Software"), to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software,
and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of
the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS
OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT
OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include <unordered_map>
#include <functional>
#include <optional>
#include <cassert>
#include <numbers>
#include <chrono>
#include <cstdio>
#include <format>
#include <memory>
#include <string>
#include <tuple>
#include <array>
#include <mutex>
#include <map>
#include <set>
#include "fft.hh"

template<bool unique, std::size_t limit>
static auto get_factors_gen(std::size_t N) // Builds a list of factors of N, up to limit. Returns the number of primes written.
{
    std::array<std::size_t,limit> target;
    static class lore /* This structure caches prime-related information. */
    {
        std::unordered_map<std::size_t, std::size_t> smallest_factors{ {0,0}, {1,1}, {2,2} };
        std::vector<signed char> is_prime{1,1,1}; // 0=not known, -1=no, 1=yes
        std::size_t              lastprime{2};    // extent of contiguous knowledge
        std::set<std::size_t>    primes{2};       // list of known primes
        std::mutex lock{};
    public:
        auto get_lock() { return std::unique_lock<std::mutex>(lock); }
        std::size_t find(std::size_t N, std::unique_lock<std::mutex>&)
        {
            if(auto i = smallest_factors.find(N); i != smallest_factors.end()) return i->second;
            std::size_t solution = N;
            for(auto p: primes)
            {
                if(N % p == 0) { solution = p; break; }
                if(p*p >= N) { break; }
            }
            for(std::size_t p = lastprime|1; p*p <= N; p+=2) if(N % p == 0) { solution = p; break; }
            if(is_prime.size() <= N) is_prime.resize(N+1);
            for(is_prime[N] = (solution==N) ? 1 : -1; lastprime+1 < is_prime.size() && is_prime[lastprime+1] != 0; ) ++lastprime;
            if(solution == N) primes.insert(solution);
            return smallest_factors.emplace(N, solution).first->second;
        }
        lore() { auto l = get_lock(); for(std::size_t a=0; a<1048576*2; ++a) { find(a, l); if(a<31) find(1<<a, l); } }
    } prime_lore;
    auto lock = prime_lore.get_lock();
    for(std::size_t orig=N,count=0,prev=0,f; ; N/=f)
        if(f=prime_lore.find(N, lock); count >= limit || f <= 1 || f == orig) { return std::pair(target,count); }
        else if(f != prev || !unique) { target[count++] = prev = f; }
}

// log2 of the product of first 16 primes ≈ 64.82. If size_t ≤ 64 bits, then size 16 array is enough.
auto& get_factors_unique = get_factors_gen<true, 16>;

// The smallest prime is 2. The largest size_t is 2^64, which has 64 factors of 2.
auto& get_factors_all    = get_factors_gen<false, 64>;

static const float pi = std::numbers::pi_v<float>;
static complex w(std::size_t k, std::size_t N) { return std::polar(1.0f, (-2 * pi) * k / N); }
static constexpr complex i(0,1);


template<typename T>
struct like_ref_array
{
    T* const ptr;
    std::size_t stride;
    T& operator[](std::size_t n) { return ptr[n*stride]; }
};

template<std::size_t N>
struct complex_array
{
    using me = complex_array<N>;
    std::array<complex,N> val;
    complex_array() = default;
    complex_array(const me&) = default;
    complex_array(const std::array<complex,N>& b) : val(b) {}
    inline me& operator=(const me& b)
    {
        for(std::size_t n=0; n<N; ++n) val[n]=b.val[n];
        return *this;
    }
    inline me& operator=(const complex& b) { for(auto& v: val) v = b; return *this; }
    template<typename U> inline me& operator=(U&& b) { val=b.val; return *this; }
    inline me operator+(const me& b) const
    {
        std::array<complex,N> res;
        for(std::size_t n=0; n<N; ++n) res[n] = val[n] + b.val[n];
        return {res};
    }
    inline me operator-(const me& b) const
    {
        std::array<complex,N> res;
        for(std::size_t n=0; n<N; ++n) res[n] = val[n] - b.val[n];
        return {res};
    }
    inline me operator*(float b) const
    {
        std::array<complex,N> res;
        for(std::size_t n=0; n<N; ++n) res[n] = val[n] * b;
        return {res};
    }
    inline me operator*(const complex& b) const
    {
        std::array<complex,N> res;
        for(std::size_t n=0; n<N; ++n) res[n] = val[n] * b;
        return {res};
    }
};
#define sab(a,b) std::sin(pi*a/b)
#define cab(a,b) std::cos(pi*a/b)
#define BEGIN_DFT_METHOD(name,is_inplace,tpl,getname) \
    tpl struct DFT_##name: public DFTbase \
    { \
        DFT_##name(long n, std::monostate={}): DFTbase(n,is_inplace){} \
        virtual std::string get_name() const override { return getname; } \
        virtual void xform_many(const complex* in, std::size_t istep,std::size_t istep2, \
                                complex* out, std::size_t ostep,std::size_t ostep2, std::size_t num) override \
        { \
            if(!num) [[unlikely]] return; \
            if(num==1)      xform_impl<false>(in,istep,istep2,out,ostep,ostep2,num); \
            else [[likely]] xform_impl<true>(in,istep,istep2,out,ostep,ostep2,num); \
        } \
        template<bool many> \
        void xform_impl(const complex* in, std::size_t istep,std::size_t istep2, \
                        complex* out, std::size_t ostep,std::size_t ostep2, std::size_t num) \
        { \
            if constexpr(!many) { assert(istep==1); assert(ostep==1); istep=ostep=1; istep2=ostep2=0; num=1; }
#define END_DFT_METHOD \
        } \
    }

#define NEXT_DFT_METHOD(name,is_inplace,tpl,nm)      END_DFT_METHOD; BEGIN_DFT_METHOD(name,is_inplace,tpl,nm)
#define NEXT_DFT_METHOD_LOOP(name,is_inplace,tpl,nm) END_DFT_METHOD; BEGIN_DFT_METHOD(name,is_inplace,tpl,nm) \
        constexpr std::size_t N=name, Z=~std::size_t{}; \
        if constexpr(!(simd&1)) goto r; \
        /*if(num==1) [[likely]] goto r;*/ \
        if(istep == 1) \
        { \
            if(ostep == 1) \
                if(istep2 == N && ostep2 == N) \
                    dosimd2<1,N,1,N>(in,istep,istep2,out,ostep,ostep2, num); \
                else \
                    dosimd2<1,Z,1,Z>(in,istep,istep2,out,ostep,ostep2, num); \
            else \
                dosimd2<1,Z,Z,Z>(in,istep,istep2,out,ostep,ostep2, num); \
        } else if(ostep == 1) \
        { \
            if(istep2 == N && ostep2 == N) \
                dosimd2<Z,N,1,N>(in,istep,istep2,out,ostep,ostep2, num); \
            else \
                dosimd2<Z,Z,1,Z>(in,istep,istep2,out,ostep,ostep2, num); \
        } \
        else if(istep2 == 1 && ostep2 == 1) [[likely]] \
            dosimd2<Z,1,Z,1>(in,istep,istep2,out,ostep,ostep2, num); \
        else r:\
            dosimd2(in,istep,istep2,out,ostep,ostep2, num); \
    } \
    template<std::size_t IS=~std::size_t{},std::size_t IS2=~std::size_t{}, \
             std::size_t OS=~std::size_t{},std::size_t OS2=~std::size_t{}> static inline \
    void dosimd2(const complex* in,std::size_t istep,std::size_t istep2, \
                 complex*      out,std::size_t ostep,std::size_t ostep2, std::size_t num) \
    { \
        if constexpr(IS  != ~std::size_t{}) istep  = IS; \
        if constexpr(IS2 != ~std::size_t{}) istep2 = IS2; \
        if constexpr(OS  != ~std::size_t{}) ostep  = OS; \
        if constexpr(OS2 != ~std::size_t{}) ostep2 = OS2; \
        std::size_t n=0; \
        if constexpr(simd>=8) \
        for(; (n+8) <= num; n+=8) { dosimd3<8,IS,IS2,OS,OS2>(in,istep,istep2, out,ostep,ostep2); in+=8*istep2; out+=8*ostep2; } \
        if constexpr(simd>=4) if(num >= 8) [[unlikely]] \
        for(; (n+4) <= num; n+=4) { dosimd3<4,IS,IS2,OS,OS2>(in,istep,istep2, out,ostep,ostep2); in+=4*istep2; out+=4*ostep2; } \
        for(; (n+1) <= num; n+=1) { dosimd3<1,IS,IS2,OS,OS2>(in,istep,istep2, out,ostep,ostep2); in+=1*istep2; out+=1*ostep2; } \
    } \
    template<std::size_t U, std::size_t IS=~std::size_t{},std::size_t IS2=~std::size_t{}, \
                            std::size_t OS=~std::size_t{},std::size_t OS2=~std::size_t{}> static inline \
    void dosimd3(const complex* in,std::size_t istep,std::size_t istep2, \
                 complex*      out,std::size_t ostep,std::size_t ostep2) \
    { \
        constexpr std::size_t N=name; \
        if constexpr(IS  != ~std::size_t{}) istep  = IS; \
        if constexpr(IS2 != ~std::size_t{}) istep2 = IS2; \
        if constexpr(OS  != ~std::size_t{}) ostep  = OS; \
        if constexpr(OS2 != ~std::size_t{}) ostep2 = OS2; \
        if constexpr(U==1) \
        { \
            like_ref_array x{in,istep}; \
            like_ref_array X{out,ostep}; \
            transform<complex>(x,X); \
        } \
        else \
        { \
            complex_array<U> x[N], X[N]; \
            for(std::size_t a=0; a<N; ++a) \
                for(std::size_t b=0; b<U; ++b) \
                    x[a].val[b] = in[a*istep + b*istep2]; \
            transform<complex_array<U>>(x, X); \
            for(std::size_t a=0; a<N; ++a) \
                for(std::size_t b=0; b<U; ++b) \
                    out[a*ostep + b*ostep2] = X[a].val[b]; \
        } \
    } \
    template<typename complex, typename z,typename Z> static inline \
    void transform(z& x, Z& X) \
    {

#define DFT_CONSTRUCT(name, params, ...) DFT_##name params: __VA_ARGS__ {
#define DFT_PROPERTIES }

#define x₀ x[0]
#define x₁ x[1]
#define x₂ x[2]
#define x₃ x[3]
#define x₄ x[4]
#define x₅ x[5]
#define x₆ x[6]
#define x₇ x[7]
#define X₀ X[0]
#define X₁ X[1]
#define X₂ X[2]
#define X₃ X[3]
#define X₄ X[4]
#define X₅ X[5]
#define X₆ X[6]
#define X₇ X[7]

// Generic fall-back method for *any* size
// Discrete Fourier Transform: Xₙ = ∑ₖ₌₀ᴺ⁻¹ (xₖ · exp(-2iπnk/N))
BEGIN_DFT_METHOD(DFT,false,, "DFT("+std::to_string(N)+')')
    std::vector<complex> tmp(num);
    for(std::size_t a=0; a<N; ++a)
    {
        for(std::size_t n=0; n<num; ++n) tmp[n] = in[0*istep + n*istep2];
        for(std::size_t b=1; b<N; ++b)
        {
            complex s = w(a*b, N);
            for(std::size_t n=0; n<num; ++n) tmp[n] += in[b*istep + n*istep2] * s;
        }
        for(std::size_t n=0; n<num; ++n) out[a*ostep + n*ostep2] = tmp[n];
    }

NEXT_DFT_METHOD_LOOP(1,true,template<int simd>,"1_"+std::to_string(simd))
    X₀ = x₀;

NEXT_DFT_METHOD_LOOP(2,true,template<int simd>,"2_"+std::to_string(simd))
    std::tie(X₀,X₁) = std::tuple(x₀+x₁, x₀-x₁);

NEXT_DFT_METHOD_LOOP(3,true,template<int simd>,"3_"+std::to_string(simd))
    const complex t₀=(x₁-x₂)*sab(1,3)*i;
    const complex u₀=(x₁+x₂),  u₁=x₀-u₀*.5f;
    std::tie(X₀,X₁,X₂)=std::tuple(x₀+u₀, u₁-t₀, u₁+t₀);

NEXT_DFT_METHOD_LOOP(4,true,template<int simd>,"4_"+std::to_string(simd))
    const complex t₀=x₀+x₂, t₁=(x₃+x₁);
    const complex u₀=x₀-x₂, u₁=(x₃-x₁)*i;
    std::tie(X₀,X₁,X₂,X₃) = std::tuple(t₀+t₁, u₀+u₁,
                                       t₀-t₁, u₀-u₁);
NEXT_DFT_METHOD_LOOP(5,true,template<int simd>,"5_"+std::to_string(simd))
    const float a=.25f, b=sab(2,5), c=sab(1,5), d=cab(1,5)-a;
    const complex t₀=x₁+x₄, t₁=x₂+x₃, t₂=(t₀-t₁)*d;
    const complex u₀=x₁-x₄, u₁=x₂-x₃, u₂=(t₀+t₁), u₃=x₀-u₂*a;
    const complex t₄=u₃+t₂, t₅=(u₀*b + u₁*c)*i;   X₀=x₀+u₂;
    const complex u₄=u₃-t₂, u₅=(u₁*b - u₀*c)*i;
    std::tie(X₁,X₂,X₄,X₃) = std::tuple(t₄-t₅, u₄+u₅,
                                       t₄+t₅, u₄-u₅);
NEXT_DFT_METHOD_LOOP(6,true,template<int simd>,"6_"+std::to_string(simd))
    const float a=.5f, b=sab(1,3);
    const complex t₀=x₀+x₃, t₁=x₄+x₁, t₂=x₂+x₅, t₃=t₀-(t₁+t₂)*a,   t₄=(t₁-t₂)*i;
    const complex u₀=x₀-x₃, u₁=x₄-x₁, u₂=x₂-x₅, u₃=u₀-(u₁+u₂)*a*i, u₄=(u₁-u₂)*i;
    std::tie(X₀,X₁,X₄,X₃,X₅,X₂) = std::tuple(t₀+t₁+t₂, u₃+u₄*b, t₃+t₄*b,
                                             u₀+u₁+u₂, u₃-u₄*b, t₃-t₄*b);
NEXT_DFT_METHOD_LOOP(8,true,template<int simd>,"8_"+std::to_string(simd))
    const float a = sab(1,4);
    const complex t₀=x₇-x₃, t₁=x₁-x₅, t₂=x₀+x₄, t₃=x₂+x₆, t₄=(t₀+t₁)*a;
    const complex u₀=x₇+x₃, u₁=x₁+x₅, u₂=x₀-x₄, u₃=x₂-x₆, u₄=(t₀-t₁)*a;
    const complex t₅=t₂+t₃, t₆=u₂+t₄, t₇=(u₀+u₁),   t₈=(u₄-u₃)*i;
    const complex u₅=t₂-t₃, u₆=u₂-t₄, u₇=(u₀-u₁)*i, u₈=(u₄+u₃)*i;
    std::tie(X₀,X₁,X₂,X₃,X₄,X₇,X₆,X₅) = std::tuple(t₅+t₇,t₆+t₈,u₅+u₇,u₆+u₈,
                                                   t₅-t₇,t₆-t₈,u₅-u₇,u₆-u₈);

// Cooley-Tukey — A fast method for radix-P sizes
NEXT_DFT_METHOD(RadixP,false, template<std::size_t P>, std::format("RadixP<{}>({})", P,N))
    std::size_t Q = N/P; // P*Q = N
    /* P transforms of Q size each */
    if(istep == istep2 && ostep*Q == ostep2)
    {
        DFTbase* DFTq = FindDFT(Q, P*istep,istep, ostep,ostep*Q, P*num);
        DFTq->xform_many(in,P*istep,istep,  out,ostep,ostep*Q, P*num);
    }
    else
    {
        DFTbase* DFTq = FindDFT(Q, P*istep,istep, ostep,ostep*Q, P);
        for(std::size_t n=0; n<num; ++n)
            DFTq->xform_many(in + n*istep2,P*istep,istep,  out + n*ostep2,ostep,ostep*Q, P);
    }
    for(std::size_t b=1; b<P; ++b)
        for(std::size_t a=1; a<Q; ++a)
            for(std::size_t n=0; n<num; ++n)
                out[ostep*(b*Q+a) + n*ostep2] *= wtable[b*Q+a];
    /* Q transforms of P size each */
    if(ostep == ostep2)
    {
        DFTbase* DFTm = FindDFT(P, ostep*Q,ostep, ostep*Q,ostep,  Q*num);
        DFTm->xform_inplace_many(out,ostep*Q,ostep, Q*num);
    }
    else
    {
        DFTbase* DFTp = FindDFT(P, ostep*Q,ostep, ostep*Q,ostep,  Q);
        for(std::size_t n=0; n<num; ++n)
            DFTp->xform_inplace_many(out + n*ostep2,ostep*Q,ostep, Q);
    }
  DFT_PROPERTIES
    std::vector<complex> wtable;
  DFT_CONSTRUCT(RadixP, (std::size_t N), DFTbase{N,false}, wtable(N))
    for(std::size_t Q=N/P, a=0; a<N; ++a)
        wtable[a] = w((a%Q)*(a/Q), N);

// Rader — A good method for prime sizes
NEXT_DFT_METHOD(Rader,true,, std::format("Rader({})", N))
    DFTbase* DFTn1 = FindDFT(N-1, 1,N-1, 1,N-1, num);

    bool needs_2sets = !DFTn1->inplace;
    std::vector<complex> buf(num + (N-1)*num*(needs_2sets?2:1));
    complex *in0 = &buf[0], *set1 = &buf[num], *set2 = &set1[needs_2sets ? (N-1)*num : 0];
    // If DFTn1 can work efficiently in-place, we can save memory and work with 1 buffer instead of 2.
    for(std::size_t n=0; n<num; ++n) in0[n] = in[0*istep + n*istep2];
    // Permutation of input
    for(std::size_t gₚ=1, k=0; k<N-1; ++k, gₚ = (gₚ*g)%N)
        for(std::size_t n=0; n<num; ++n)
            set1[k+n*(N-1)] = in[gₚ*istep + n*istep2];
    // DFT
    DFTn1->xform_many(&set1[0],1,N-1, &set2[0],1,N-1, num);
    // Set output DC component:
    for(std::size_t n=0; n<num; ++n) out[0*ostep + n*ostep2] = in0[n] + set2[n*(N-1)];
    // Multiply by ω:
    for(std::size_t n=0; n<num; ++n)
        for(std::size_t k=0; k<N-1; ++k)
            set2[k+(N-1)*n] = std::conj( set2[k+(N-1)*n] * ω[k] + (k==0 ? in0[n] : 0) );
    // IDFT
    DFTn1->xform_many(&set2[0],1,N-1, &set1[0],1,N-1, num);
    // Inverse permutation
    for(std::size_t gₚ=1, k=0; k<N-1; ++k, gₚ = (gₚ*gᵢₙᵥ)%N)
        for(std::size_t n=0; n<num; ++n)
            out[gₚ*ostep + n*ostep2] = std::conj(set1[k+(N-1)*n]);
  DFT_PROPERTIES
    std::size_t g, gᵢₙᵥ;
    std::vector<complex> ω;
    /* Compute n^m mod p, where m >= 0 and p > 0. */
    static std::size_t powermod(std::size_t n, std::size_t m, std::size_t p)
    {
        if(m == 0) return 1;
        if(m % 2 == 0) { n = powermod(n, m/2, p); return (n*n) % p; }
        return (n * powermod(n, m-1, p)) % p;
    }
  DFT_CONSTRUCT(Rader, (std::size_t N), DFTbase{N,true}, g{2}, ω(N-1))
    auto [factors, count] = get_factors_unique(N-1);
    for(std::size_t m=0; m<count; )
        if(powermod(g, (N-1) / factors[m++], N) == 1)
            { ++g; m = 0; }
    gᵢₙᵥ = powermod(g,N-2,N);
    for(std::size_t gₚ=1, i=0; i < N-1; ++i, gₚ = (gₚ*gᵢₙᵥ)%N)
        { ω[i] = w(gₚ, N); }
    FindDFT(N-1, 1,1, 1,1, 1)->xform_inplace(&ω[0], 1);
    for(std::size_t i=0; i<N-1; ++i) ω[i] /= float(N-1);

// Bluestein — An ok method (sometimes fast too) for all sizes
NEXT_DFT_METHOD(Bluestein,true,, std::format("Bluestein({})(nb={})", N, nb))
    DFTbase* DFTnb = FindDFT(nb,1,nb, 1,nb, num);
    bool needs_2sets = !DFTnb->inplace;
    std::vector<complex> buf(nb*num*(needs_2sets ? 2 : 1));
    complex *set1 = &buf[0], *set2 = &set1[needs_2sets ? nb*num : 0];
    // If DFTnb can work efficiently in-place, we can save memory and work with 1 buffer instead of 2.
    for(std::size_t n=0; n<num; ++n)
    for(std::size_t k=0; k<N; ++k) { set1[k+n*nb] = in[k*istep+n*istep2] * ω₀[k]; }
    DFTnb->xform_many(&set1[0],1,nb, &set2[0],1,nb, num); // DFT 
    for(std::size_t n=0; n<num; ++n)
    for(std::size_t j=0; j<nb; ++j) { set2[j+n*nb] = std::conj(set2[j+n*nb]) * ω₁[j]; }
    DFTnb->xform_many(&set2[0],1,nb, &set1[0],1,nb, num); // IDFT
    for(std::size_t n=0; n<num; ++n)
    for(std::size_t k=0; k<N; ++k) { out[k*ostep+n*ostep2] = std::conj(set1[k+n*nb]) * ω₀[k]; }
  DFT_PROPERTIES
    std::size_t nb;
    std::vector<complex> ω₀, ω₁;
  DFT_CONSTRUCT(Bluestein, (std::size_t N, std::size_t Nb), DFTbase{N,true}, nb{Nb}, ω₀(N), ω₁(nb))
    assert(nb >= 2*N-1);
    for(std::size_t k=0; k<N; ++k) { ω₀[k] = w(k*k, 2*N); }
    for(std::size_t k=0; k<N; ++k) { ω₁[k] = ω₀[k] / float(nb); }
    for(std::size_t k=1; k<N; ++k) { ω₁[nb-k] = ω₁[k]; }
    FindDFT(nb, 1,1, 1,1, 1)->xform_inplace(&ω₁[0],1);
END_DFT_METHOD;


DFTbase* FindDFT(std::size_t N, std::size_t istep,std::size_t istep2,
                                std::size_t ostep,std::size_t ostep2,
                                std::size_t num)
{
    // Note: This function doesn’t use istep,ostep,istep2,ostep2,num parameters
    // for anything, but if one wants to do change the function so that it
    // chooses the method by doing performance comparisons, these would be
    // valuable.
    if(num == 1) istep2=ostep2=1;

    // Use memoization for FFT methods.
    static std::pair<std::mutex,std::map<std::size_t, std::unique_ptr<DFTbase>>> data;
    auto& [lock,table] = data;
    if(auto i = table.find(N); i != table.end()) return &*i->second;
    // Perform a second search under exclusive access
    std::unique_lock lk(lock);
    if(auto i = table.find(N); i != table.end()) return &*i->second;
    lk.unlock();
    //std::fprintf(stderr, "FindDFT(%zu)\n", N);
    DFTbase* res = [=]->DFTbase*{
        switch(N)
        {
            case 1: return new DFT_1<0>(N);
            case 2: return new DFT_2<0>(N);
            case 3: return new DFT_3<0>(N);
            case 4: return new DFT_4<0>(N);
            case 5: return new DFT_5<0>(N);
            case 6: return new DFT_6<0>(N);
            case 8: return new DFT_8<0>(N);
            case 7: return new DFT_DFT(N);
        }
        auto [factors, count] = get_factors_all(N);
        if(count >= 2)
        {
            //if(N>16 && N%16==0) return new DFT_RadixP<16>(N);
            if(N>8 && N%8==0) return new DFT_RadixP<8>(N);
            //if(N>6 && N%6==0) return new DFT_RadixP<6>(N);
            if(N>4 && N%4==0) return new DFT_RadixP<4>(N);
            #define primes(o) o(2)o(3)o(5)o(7)o(11)o(13)o(17)o(19)o(23)o(29)o(31)o(37)o(41)o(43)o(47)o(53)o(59)o(61)o(67)o(71)o(73)o(79)o(83)o(89)o(97)o(101)o(103)o(107)o(109)o(113)o(127)o(131)o(137)o(139)o(149)o(151)o(157)o(163)o(167)o(173)o(179)o(181)o(191)o(193)o(197)o(199)o(211)o(223)o(227)o(229)o(233)o(239)o(241)o(251)o(257)o(263)o(269)o(271)o(277)o(281)o(283)o(293)o(307)o(311)o(313)o(317)o(331)o(337)o(347)o(349)o(353)o(359)o(367)o(373)o(379)o(383)o(389)o(397)o(401)o(409)o(419)o(421)o(431)o(433)o(439)o(443)o(449)o(457)o(461)o(463)o(467)o(479)o(487)o(491)o(499)o(503)o(509)o(521)o(523)o(541)o(547)o(557)o(563)o(569)o(571)o(577)o(587)o(593)o(599)o(601)o(607)o(613)o(617)o(619)o(631)o(641)o(643)o(647)o(653)o(659)o(661)o(673)o(677)o(683)o(691)o(701)o(709)o(719)o(727)o(733)o(739)o(743)o(751)o(757)o(761)o(769)o(773)o(787)o(797)o(809)o(811)o(821)o(823)o(827)o(829)o(839)o(853)o(857)o(859)o(863)o(877)o(881)o(883)o(887)o(907)o(911)o(919)o(929)o(937)o(941)o(947)o(953)o(967)o(971)o(977)o(983)o(991)o(997)o(1009)o(1013)o(1019)o(1021)o(1031)o(1033)o(1039)o(1049)o(1051)o(1061)o(1063)o(1069)o(1087)o(1091)o(1093)o(1097)o(1103)o(1109)o(1117)o(1123)o(1129)o(1151)o(1153)o(1163)o(1171)o(1181)o(1187)o(1193)o(1201)o(1213)o(1217)o(1223)
            switch(factors[0])
            {
                #define o(n) case n: return new DFT_RadixP<n>(N);
                primes(o)
                #undef o
            }
        }
        // Round N up to next power of two
        std::size_t nb = 2*N-1;
        for(unsigned a=1; a< 8*sizeof(nb); a<<=1) { nb = 1 + ((nb-1) | ((nb-1) >> a)); }
        // Prime, and unwieldly for Bluestein? Use Rader.
        if(!count /*&& nb > N*5/2*/) return new DFT_Rader(N);
        // Else fallback to Bluestein.
        return new DFT_Bluestein(N, nb);
    }();
    lk.lock();
    return &*(table[N] = std::unique_ptr<DFTbase>(res));
}