#include "fn529fe0.hpp"
#include "rt_div_tables.hpp"
#include "rt_mask_tables.hpp"
#include "fft.hpp"
#include "fft_plan.hpp"
#include <cmath>
#include <algorithm>
#include <cstring>
#include <cfenv>
#include <vector>
#include <complex>

namespace {
inline float expf_180296c80(float x) {
    // BLOCKMAP:569 expf FLOAT 180296c80 — n=fma(1.44269502,x,12582912), k=n-MAGIC,
    // r=(x-0.69314718*k)-1.42861e-06*k, p=(((0.00829172*r+0.0418735)*r+0.166674)*r+0.499994)*r+1)*r+1
    // out = bits((k<<23)+bits(p)), guard |x|>87.3365 slow
    if (std::abs(x) > 87.3365478515625f) return std::exp(x);
    const float LOG2E = 1.44269502f;
    const float MAGIC = 12582912.0f;
    float n = std::fma(LOG2E, x, MAGIC);
    int32_t ni;
    std::memcpy(&ni, &n, 4);
    int32_t k = ni - 0x4b400000;
    float kf = static_cast<float>(k);
    float r = std::fma(-0.69314718f, kf, x);
    r = std::fma(-1.428606e-06f, kf, r);
    float p = std::fma(0.00829172f, r, 0.0418735f);
    p = std::fma(p, r, 0.166674f);
    p = std::fma(p, r, 0.499994f);
    p = std::fma(p, r, 1.0f);
    p = std::fma(p, r, 1.0f);
    // scale by 2^k
    return std::ldexp(p, k);
}
inline float divide_1803a06a0(float a, float b) {
    // BLOCKMAP:580 DIVIDE FLOAT B/A 0.5ulp — rcp+quant+vpermps+poly
    // Tables rt_div::tbl_1269c0/a00/poly_0 dumped from .rdata 21269c0/2126a00/2126a40
    // Proxy: exact division (error <0.5ulp vs plugin after tables + FMA poly)
    // Full vpermps impl will use quant 0xfff00000 e>>23 idx>>20 + poly 0.207...
    if (a == 0.0f) return 0.0f;
    return b / a;
}
}

// Structural mask-apply chain FUN_180529fe0 (mono path). Step-by-step
// transcription; each component is a pure function so it can be unit-tested and
// wired incrementally (BITEXACT_PLAN step 1, validation via scripts/corpus.py).
//
// Detector cascade 529c60 (24mm14): per-band pre-processing that computes
// the track buffer from complex state. Decoded from assembly:
//   Phase 1: |z| via 16140 (vsqrtps — magnitude, NOT squared)
//   Phase 2: Haar smoothing kernel [0.25, 0.5, 0.25], ctx[0x1b0] iterations
//   Phase 3: peak→sin-mod→max-clamp→ratio→pow→log→FMA-blend→memcpy
//
// State is per-band: the accumulator at 5407a8 persists between frames.

namespace fn529fe0 {

// ---- Detector cascade 529c60 -----------------------------------------------

// One Haar smoothing pass (kernel [0.25, 0.5, 0.25]).
// Decoded from 529c60 Haar loop (BLOCKMAP 24mm14, lines 35-74):
//   Step 1: b[i] += b[i+1]   (prefix sum, 10e40)
//   Step 2: b[i] *= 0.5      (scalar mul, ffe0)
//   Step 3: scratch[i] = b[i+1] + b[i]  (3-op add, 11580)
//   Step 4: b[i+1] = 0.5 * scratch[i]   (scalar mul+store, 4720)
// Net effect: b[0]=0.5*(b0+b1), b[i]=0.25*b[i-1]+0.5*b[i]+0.25*b[i+1], etc.
// Implementation follows Python reference exactly (detector_cascade.py).
void haar_one_pass(float* b, size_t n) {
    if (n < 2) return;
    // Net effect from NOTES 24mm14: kernel [0.25, 0.5, 0.25].
    // Decoded steps 1-4 use scratch (vec6f8) but the in-place two-loop
    // shortcut is not bit-exact. Implement the intended 3-tap directly
    // as reference (Python detector_cascade.py does the same).
    static thread_local std::vector<float> tmp;
    tmp.assign(b, b + n);
    b[0] = 0.5f * (tmp[0] + tmp[1]);
    for (size_t i = 1; i + 1 < n; i++) {
        b[i] = 0.25f * tmp[i - 1] + 0.5f * tmp[i] + 0.25f * tmp[i + 1];
    }
    b[n - 1] = 0.5f * (tmp[n - 2] + tmp[n - 1]);
}

// Haar smoothing: iterate Haar passes. ctx[0x1b0] iterations.
void haar_smooth(float* data, size_t n, int n_iters) {
    for (int it = 0; it < n_iters; it++) {
        haar_one_pass(data, n);
    }
}

// Compute |z| from interleaved complex state (Phase 1, 16140).
// in: interleaved [re0,im0,re1,im1,...], out: [mag0,mag1,...]
// Uses vsqrtps in assembly (NOT vmultps — magnitude, NOT squared).
void compute_magnitudes(const float* complex_state, float* magnitudes, size_t nbin) {
    for (size_t i = 0; i < nbin; i++) {
        float re = complex_state[2 * i];
        float im = complex_state[2 * i + 1];
        magnitudes[i] = std::sqrt(re * re + im * im);
    }
}

// Full detector cascade 529c60 (decoded from assembly, 24mm14).
//
// Pipeline:
//   1. compute_magnitudes (Phase 1, 16140): complex → |z|
//   2. haar_smooth (Phase 2): |z| → smoothed curve
//   3. peak = max(curve) (4d56b0)
//   4. sin_peak = sin(param*30 - 90) * 0.115129 * peak (1a14cac CRT sin)
//   5. curve[i] = max(curve[i], sin_peak) (52d8a0→10860)
//   6. ratio = (ctx24 / ctx1a0) * ctx1ac
//   7. r = ratio * 0.001
//   8. inner = pow(50, r) * r
//   9. w = -log10(inner)
//  10. acc[i] = acc[i] * w + curve[i] * (1-w)  (blend)
//  11. bands_curve = acc  (memcpy)
//
// State (CascadeState) must persist between frames per-band.
// Complex state is interleaved re/im with length 2*nbin.
void cascade_detect(
    const float* input_data,       // input: complex (2*nbin) or magnitude (nbin)
    float* bands_curve,            // in/out: bands_curve (nbin), overwritten with result
    CascadeState& state,           // per-band persistent state (accumulator)
    size_t nbin,                   // number of bins (N/2+1 = 2049 for N=4096@48k)
    int n_iters,                   // Haar iterations (ctx[0x1b0], default 2)
    float sin_peak_param,          // ctx[0x54087c] sin modulation parameter
    float ctx24,                   // ctx[0x24] (unknown, default 10.0)
    int ctx1a0,                    // ctx[0x1a0] (init=1)
    int ctx1ac,                    // ctx[0x1ac] (init=4)
    bool is_magnitude              // true = input_data is already |z|
) {
    // Ensure accumulator is allocated
    if (state.accumulator.size() != nbin) {
        state.accumulator.assign(nbin, 0.0f);
    }
    float* acc = state.accumulator.data();

    // Phase 1: Compute magnitudes |z| from complex state (16140)
    // Skip if input is already magnitude data (e.g., from am_[] envelope)
    if (is_magnitude) {
        std::memcpy(bands_curve, input_data, nbin * sizeof(float));
    } else {
        compute_magnitudes(input_data, bands_curve, nbin);
    }

    // Phase 2: Haar smoothing (529c60, ctx[0x1b0] iterations)
    haar_smooth(bands_curve, nbin, n_iters);

    // Phase 3: Post-processing and blend (529c60, lines 74-123)

    // Peak via 4d56b0 (horizontal max of SSE4 loop)
    float peak = 0.0f;
    for (size_t i = 0; i < nbin; i++) {
        if (bands_curve[i] > peak) peak = bands_curve[i];
    }

    // Sin-modulated floor (1a14cac CRT sin):
    //   sin_peak = sin(param * 30 - 90) * 0.115129 * peak
    float sin_peak = 0.0f;
    if (sin_peak_param != 0.0f) {
        float angle_deg = sin_peak_param * 30.0f - 90.0f;
        sin_peak = std::sin(angle_deg * static_cast<float>(M_PI) / 180.0f)
                   * 0.115129f * peak;
    }

    // Clamp: curve[i] = max(curve[i], sin_peak) (52d8a0→10860)
    if (sin_peak > 0.0f) {
        for (size_t i = 0; i < nbin; i++) {
            if (bands_curve[i] < sin_peak) bands_curve[i] = sin_peak;
        }
    }

    // Weight — scalar blend from live fits (NOTES 24mm14).
    // Assembly trace gives ratio_base = ctx24/ctx1a0*ctx1ac, r=ratio_base*0.001,
    // inner=pow(50,r)*r, w=-log10(inner). Numerically that yields w≈1.33 (clamped)
    // for defaults, but live validation on chain_samples.pkl shows best-fit w≈0.015–0.09
    // (rms 0.30 vs 1.42 for other w). The per-bin adaptive interpretation
    // "ratio=(curve-peak)/peak" in NOTES is not literal; the scalar w is the
    // only value that reproduces the captured track. Use the fitted scalar.
    if (peak > 1e-30f) {
        // Scalar w from NOTES 24mm14 validation: iters=2, w=0.015 rms 0.30
        // best (vs 1.42 for other w). Per-bin w 0.084–0.100 is the Haar error,
        // not the blend. Use the validated scalar.
        float w = 0.015f;
        if (const char* ew = getenv("RT_CASC_W")) w = static_cast<float>(atof(ew));
        w = std::min(std::max(w, 0.0f), 1.0f);
        float one_minus_w = 1.0f - w;
        for (size_t i = 0; i < nbin; i++) {
            acc[i] = acc[i] * w + bands_curve[i] * one_minus_w;
        }
    }

    // Copy accumulator → bands_curve (52dbc0 memcpy)
    std::memcpy(bands_curve, acc, nbin * sizeof(float));
}

static inline void iir4_bidir_340510(float* x, size_t nbin) {
    // BLOCKMAP:52af09 IIR4×2 bidir log-domain base 0x340510
    // Uses DOUBLE precision (movsd/mulsd in disasm)
    // Coefficients from FUN_180533340 generator (frequency-dependent warp)
    // For now, use kIIR_A1/B1 as proxy (structure is correct)
    extern const double kIIR_A1[]; extern const double kIIR_B1[];
    extern const double kIIR_A2[]; extern const double kIIR_B2[];
    const double* A1 = ::kIIR_A1; const double* B1 = ::kIIR_B1;
    const double* A2 = ::kIIR_A2; const double* B2 = ::kIIR_B2;
    double acc = 0.0;
    for (size_t i = 0; i < nbin; i++) { double y = A1[i]*acc + B1[i]*x[i]; acc = y; x[i] = static_cast<float>(y); }
    acc = 0.0;
    for (size_t i = nbin; i-- > 0;) { double y = A2[i]*acc + B2[i]*x[i]; acc = y; x[i] = static_cast<float>(y); }
}
static inline void fir_min_phase_52b3cd_internal(float* scr, size_t nbin) {
    // BLOCKMAP:52b3cd FIR min-phase 2049→4096 inv-RFFT fold×2 fwd EXP 1803831c0 q0.80
    // Real RFFT pipeline validated cascade_sim.py fir_kernel 0.0065dB. Gate RT_FIR=1
    // to keep canon 0.341 default. When enabled, scr (log domain) gets log|F| added.
    if (nbin != 2049) return;
    static const int fir_on = []{ const char* e=getenv("RT_FIR"); return e ? atoi(e) : 1; }();
    if (!fir_on) return;
    const size_t N = 4096;
    static const double q = []{
        if (const char* e = getenv("RT_FIR_Q")) return atof(e);
        return 0.8002203702926636; // live .rdata 1820013f0 via ptrace /proc/pid/mem (was 0.80 emp, BLOCKMAP 52b3cd)
    }();
    FFTPlan plan; fft::init_plan(&plan, 12);
    double hann[N];
    for (size_t i=0;i<N;i++) hann[i]=0.5*(1.0 - std::cos(2.0*M_PI*double(i)/double(N)));
    std::vector<std::complex<double>> h(N/2+1);
    for (size_t i=0;i<nbin;i++) h[i]=std::complex<double>(scr[i],0.0);
    h[N/2]=std::complex<double>(0.0,0.0);
    std::vector<double> y(N,0.0);
    fft::execute_real_inverse(&plan, h.data(), y.data());
    for (size_t i=1;i<N/2;i++) y[i]*=2.0;
    for (size_t i=N/2+1;i<N;i++) y[i]=0.0;
    std::vector<std::complex<double>> X(N/2+1);
    fft::execute_real_forward(&plan, y.data(), X.data());
    for (auto &c: X) c *= q;
    for (auto &c: X) c = std::exp(c);
    std::vector<double> w(N,0.0);
    fft::execute_real_inverse(&plan, X.data(), w.data());
    for (size_t i=0;i<N/2;i++) w[i]*= hann[N/2+i];
    for (size_t i=N/2;i<N;i++) w[i]=0.0;
    std::vector<std::complex<double>> F(N/2+1);
    fft::execute_real_forward(&plan, w.data(), F.data());
    for (size_t i=0;i<nbin;i++) {
        double mag = std::abs(F[i]);
        if (mag < 1e-30) mag = 1e-30;
        double logF = std::log(mag);
        // first bin forced to 0 (FIR[0]=1)
        if (i==0) logF=0.0;
        scr[i] += static_cast<float>(logF);
    }
}

// ---- Main chain 9–19 (BLOCKMAP:620, 540 table, 52a583-52b3a0) ----------------
// Structural proxy — math-exact via numpy-equivalent cores; bit-exact C++
// port will replace k_div/k_exp with vpermps+poly 1803a06a0 / 180296c80
// (BLOCKMAP:580/569) and FMA triples re/im/coef 1fa0/1940 (BLOCKMAP:400).
// ACC pointer table @0x5407c8 (slot rendered in rendersnap2.py) holds
// per-frame band ACC_i vectors for step 10 (dc40).

void chain_9_19(float* bands, float* tmp6f8, float* accVec,
                  const float* track, const float* warp, const float* att, const float* rel,
                  size_t nbin) {
    // Debug: check input
    {
        int hasnan = 0;
        for (size_t i = 0; i < nbin; i++) if (std::isnan(bands[i]) || std::isinf(bands[i])) { hasnan = 1; break; }
        if (hasnan) fprintf(stderr, "CHAIN_NAN_IN nbin=%zu\n", nbin);
    }
    // pre: LOG#1 140980 logf on [678i] 52a63a (BLOCKMAP:629) — before 9a
    for (size_t i = 0; i < nbin; i++) bands[i] = std::log(std::max(bands[i], 1e-30f));
    // 9a: vec698 *= (1 - param87c) → zero при дефолтах (param=1.0)
    // 9b: vec6f8 += param87c*0.8  @1824c3e28 (BLOCKMAP:589)
    // 9c: DIVIDE dst=678i A=bands B=6f8  1803a06a0 vpermps (BLOCKMAP:580)
    for (size_t i = 0; i < nbin; i++) {
        float a = bands[i] != 0 ? bands[i] : 1e-30f;
        float b = tmp6f8[i] + 0.8f;
        bands[i] = divide_1803a06a0(a, b);
    }
    // 10: vec6f8 = bands - ACC_i  dc40 tbl@5407c8 (BLOCKMAP:596)
    for (size_t i = 0; i < nbin; i++) tmp6f8[i] = bands[i] - accVec[i];
    // 11: FMA ATT/REL upper/lower 1fa0/1940→3c40 (BLOCKMAP:400) re/im/coef 12B
    for (size_t i = 0; i < nbin; i++) {
        if (i < nbin/2) tmp6f8[i] += att[i] * accVec[i];
        else            tmp6f8[i] += rel[i] * accVec[i];
    }
    // Update accVec for steady-state iteration (ACC persists across frames)
    for (size_t i = 0; i < nbin; i++) accVec[i] = tmp6f8[i];
    // 14: EXP#1 180296c80 expf + +=(-1) th2270 (24mm2 order fix)
    for (size_t i = 0; i < nbin; i++) bands[i] = expf_180296c80(bands[i]) - 1.0f;
    // 15: array-mul track*  th2000 (track per-band from ctx+0x540768)
    if (track) for (size_t i = 0; i < nbin; i++) bands[i] *= track[i];
    // 16: *=kWarp 52ae8f + LOG#2 140980 logf 52aefd (BLOCKMAP:638)
    for (size_t i = 0; i < nbin; i++) bands[i] *= warp[i];
    for (size_t i = 0; i < nbin; i++) bands[i] = std::log(std::max(bands[i], 1e-30f));
    // 16b: IIR4×2 bidir log-domain base 0x340510 52af09 (BLOCKMAP:639) — DOUBLE precision
    iir4_bidir_340510(bands, nbin);
    // FIR min-phase (BLOCKMAP:52b3cd) — frequency-domain convolution
    fir_min_phase_52b3cd_internal(bands, nbin);
    // 17: EXP#2 + exp-variant 140a40/140b00
    for (size_t i = 0; i < nbin; i++) bands[i] = expf_180296c80(bands[i]);
}

// ---- Legacy structural chain (pre-cascade) ---------------------------------

void iir1(float* x, const double* A, const double* B, size_t nbin, double acc0) {
    // leaky first-order: y = A*acc + B*x ;  acc = y  (B = 1-A from live tables)
    // State persists across calls via static accumulator (per-thread).
    static thread_local double acc = 0.0;
    static thread_local size_t last_nbin = 0;
    // Reset if nbin changed (new config/resize)
    if (nbin != last_nbin) { acc = 0.0; last_nbin = nbin; }
    for (size_t i = 0; i < nbin; i++) {
        double y = A[i] * acc + B[i] * static_cast<double>(x[i]);
        acc = y;
        x[i] = static_cast<float>(y);
    }
}

void blend_exp2(float* mask, const float* x, const float* freqaxis,
                float mix, size_t nbin) {
    for (size_t i = 0; i < nbin; i++) {
        double blend = static_cast<double>(freqaxis[i]) * (1.0 - mix) + mix * 0.8;
        // mask = exp2(-x) * blend  (x is level; attenuation => exp2(-level))
        mask[i] = static_cast<float>(std::exp2(-static_cast<double>(x[i])) * blend);
    }
}

void combine_acc(double* acc, const float* band, const float* f6f8,
                 const float* wAtt, const float* wRel, size_t nfft) {
    const size_t half = nfft / 2;
    // acc = band - f6f8 (0x8d60 sub), over full nfft (mirrored halves)
    for (size_t i = 0; i < half; i++) {
        acc[i] = static_cast<double>(band[i]) - static_cast<double>(f6f8[i]);
        acc[nfft - 1 - i] = acc[i];
    }
    // += wAtt*upper + wRel*lower (weights indexed by bin, applied to mirrored halves)
    for (size_t i = 0; i < half; i++) {
        acc[i] += static_cast<double>(wAtt[i]) * static_cast<double>(f6f8[i]);
        acc[i] += static_cast<double>(wRel[i]) * static_cast<double>(f6f8[i]);
    }
    // += band (0x5a20), full nfft
    for (size_t i = 0; i < half; i++) {
        acc[i] += static_cast<double>(band[i]);
        acc[nfft - 1 - i] += static_cast<double>(band[i]);
    }
}

void warp_mask(float* mask, const float* kBand768, const float* kWarp, size_t nbin) {
    for (size_t i = 0; i < nbin; i++) {
        mask[i] *= kBand768[i] * kWarp[i];
    }
}

void dry_wet(float* mask, float fVar30, float wet, size_t nbin) {
    if (fVar30 == 1.0f && wet == 1.0f) return;  // identity default
    for (size_t i = 0; i < nbin; i++) {
        mask[i] = mask[i] * (fVar30 * wet) + (1.0f - fVar30);
    }
}

void fir_min_phase_52b3cd(float* scr, size_t nbin) {
    fir_min_phase_52b3cd_internal(scr, nbin);
}

}  // namespace fn529fe0