radix4.hpp

template <
        bool Inv,
        std::ptrdiff_t N_out
    >
    inline void
    butterfly_radix4(gsl::span<cpx_t, N_out> output)
    const noexcept {
        auto const m = N_out/4;
        auto const stride = (N_out > 0) ? N / N_out : 0;

        std::swap(output[1*m], output[2*m]);
        std::tie(output[0*m], output[1*m]) = math::scissors(output[0*m], output[1*m]);
        std::tie(output[2*m], output[3*m]) = math::scissors(output[2*m], output[3*m]);
        std::tie(output[0*m], output[2*m]) = math::scissors(output[0*m], output[2*m]);
        output[3*m] = flip<Inv>(output[3*m]);
        std::tie(output[1*m], output[3*m]) = math::scissors(output[1*m], output[3*m]);

        for (auto i = 1; i < m; ++i) {
            std::swap(output[1*m], output[2*m]);
            output[1*m] = math::multiply_fast(output[1*m], twiddle<Inv>(2*i * stride));
            output[2*m] = math::multiply_fast(output[2*m], twiddle<Inv>(1*i * stride));
            output[3*m] = math::multiply_fast(output[3*m], twiddle<Inv>(3*i * stride));
            std::tie(output[0*m], output[1*m]) = math::scissors(output[0*m], output[1*m]);
            std::tie(output[2*m], output[3*m]) = math::scissors(output[2*m], output[3*m]);
            std::tie(output[0*m], output[2*m]) = math::scissors(output[0*m], output[2*m]);
            output[3*m] = flip<Inv>(output[3*m]);
            std::tie(output[1*m], output[3*m]) = math::scissors(output[1*m], output[3*m]);
        }
    }