#pragma once

// Bit-exact transcription of the sincospi "rotor" kernel 0x180184900 (fast path),
// per-element scalar form, mirroring the exact VEX/FMA op order of big_184900.dis
// (main loop 0x180184a80-0x1801850d4 / masked tail 0x180185160-0x18018530c).
//
// Per input float x the kernel emits two floats:
//     bufA[i] = sin(x)   (buf1=rdx/arg at kernel entry)
//     bufB[i] = cos(x)   (buf2=r8 /arg)
// (verified: tail loop stores sin->[r11+rbx]=bufA via ymm5, cos->[r11+r12]=bufB via ymm2;
//  main loop stores the 4 lane-group results ymm{7,6,4,1}->bufA, ymm{8,5,3,2}->bufB).
//
// The 8-lane AVX body interleaves four independent groups; all ops are per-lane
// (no cross-lane shuffle), so the scalar loop below reproduces each element exactly.
//
// The slow path (0x180186086, entered when any lane |x| >= 10000 or non-finite) and the
// inf/nan guard helper 0x180189d60 are NOT reproduced: for audio |x| <= 1 they never run.
// Contract: inputs must satisfy |x| < 10000 (checked in rotor_sincos via nan guard).

#include <cstdint>
#include <cstring>
#include <cmath>

namespace detkernel {

// exact fused multiply-add with a single rounding (no re-contraction by the compiler)
static inline float fma_f(float a, float b, float c) {
    return __builtin_fmaf(a, b, c);
}

static inline uint32_t bit_cast_u32(float f) {
    uint32_t u;
    std::memcpy(&u, &f, sizeof(u));
    return u;
}

static inline float bit_cast_f(uint32_t u) {
    float f;
    std::memcpy(&f, &u, sizeof(f));
    return f;
}

// magic-rounding trick constants (from tables 0x181d1b3x0 / 0x181d1b0c0-0x181d1b300)
static constexpr uint32_t C_ABS_MASK   = 0x7fffffffu;   // 0x181d1ad40
static constexpr float    C_INV_PI     = 0.318309873f;  // 0x181d1b340  4b400000? no: 3ea2f983
static constexpr float    C_MAGIC      = 12582912.0f;   // 0x181d1b380  0x4b400000
static constexpr float    C_PI_HI      = 3.141592741f;  // 0x181d1b0c0  0x40490fdb
static constexpr float    C_PI_LO      = -8.742277658e-08f; // 0x181d1b100 0xb3bbbd2e
static constexpr float    C_RED_CORR   = -3.430249024e-15f; // 0x181d1b140 0xa7772ced (C0)
static constexpr float    C_SIN_1      = 0.00833306462f;    // 0x181d1b280 0x3c088768 (C1)
static constexpr float    C_SIN_2      = -1.980916067e-04f; // 0x181d1b2c0 0xb94fb6cf (C2)
static constexpr float    C_SIN_3      = 2.604164590e-06f;  // 0x181d1b300 0x362ec335 (C3)
static constexpr float    C_SIN_N6     = -0.166666612f;     // 0x181d1b180 0xbe2aaaa7
static constexpr uint32_t C_HALF_BITS  = 0x3f000000u;       // 0.5f  0x181d1b400
static constexpr float    C_ONE        = 1.0f;              // 0x181d1b440
static constexpr float    C_NOISE      = 10000.0f;          // 0x181d1ad80 (slow-path threshold)
static constexpr float    C_MASK_DEF   = 0.75f;             // 0x181d1c280 (filled default lanes)

// rotor_sincos(x, &sin, &cos) per 0x180184900 (elementwise, exact op order)
static inline void rotor_sincos(float x, float& s, float& c) {
    uint32_t ux;
    std::memcpy(&ux, &x, sizeof(ux));

    // vmovups/ymm load + vandps abs mask + vfmadd231ps
    uint32_t ua = ux & C_ABS_MASK;
    float a;
    std::memcpy(&a, &ua, sizeof(a));

    float t = fma_f(a, C_INV_PI, C_MAGIC);   // t = magic + a*(1/pi)
    float k = t - C_MAGIC;                   // vsubps k = t - magic
    uint32_t tk;
    std::memcpy(&tk, &t, sizeof(tk));
    uint32_t sbt = tk & 1u;                  // vpslld t,31 -> bit0 (parity/round)

    float p = fma_f(-C_PI_HI, k, a);         // vfnmadd231 p = a - pi_hi*k
    p       = fma_f(-C_PI_LO, k, p);         // p -= pi_lo*k
    uint32_t up;
    std::memcpy(&up, &p, sizeof(up));
    uint32_t sbp = up & 0x80000000u;         // vandps -0.0 -> sign bit of p

    uint32_t halfbits = sbp ^ C_HALF_BITS;   // vxorps (-0.0&p) ^ 0.5  -> +-0.5
    float h;
    std::memcpy(&h, &halfbits, sizeof(h));
    float q = k + h;                         // vaddps q = k +- 0.5

    // vfnmadd213ps dest,src1,src2 = -(dest*src1) + src2
    float cc = fma_f(C_RED_CORR, -k, p);     // cc = p - C0*k          (sin arg core)
    float rq = fma_f(-C_PI_HI, q, a);        // rq = a - pi_hi*q
    rq        = fma_f(-C_PI_LO, q, rq);
    float cq = fma_f(C_RED_CORR, -q, rq);    // cq = rq - C0*q          (cos arg core)

    uint32_t sbth = sbt ? 0x80000000u : 0u;  // vpslld fully left -> 0x80000000/0
    uint32_t us1 = bit_cast_u32(cc) ^ sbth;  // arg_sin = cc ^ sbt
    // vxorps(-0.0,sbp)=0x80000000^sbp, then ^sbt -> arg_cos = cq ^ sbt ^ sbp ^ 0x80000000
    uint32_t us2 = bit_cast_u32(cq) ^ sbth ^ sbp ^ 0x80000000u;

    float sarg = bit_cast_f(us1);
    float qarg = bit_cast_f(us2);

    float ss = sarg * sarg;                  // vmulps
    float qs = qarg * qarg;

    float as = fma_f(C_SIN_3, ss, C_SIN_2);  // Horner (vfmadd231 then vfmadd213 chain)
    float aq = fma_f(C_SIN_3, qs, C_SIN_2);
    as = fma_f(as, ss, C_SIN_1);
    aq = fma_f(aq, qs, C_SIN_1);
    as = fma_f(as, ss, C_SIN_N6);
    aq = fma_f(aq, qs, C_SIN_N6);

    float os = ss * as;                      // vmulps
    float oq = qs * aq;
    float sins = fma_f(os, sarg, sarg);      // vfmadd213 sin = os*sarg + sarg
    float sinq = fma_f(oq, qarg, qarg);

    uint32_t usrc = ux & 0x80000000u;        // vandnps ~abs & src -> sign bit of src
    s = bit_cast_f(bit_cast_u32(sins) ^ usrc); // vxorps sin ^ sign(src)

    bool zero = (ux == 0u) || (ux == 0x80000000u); // vcmpeqps src == signbit(src)
    c = zero ? C_ONE : sinq;                  // vblendvps -> 1.0 for zero lanes
}

// 8-lane rotor; mirrors the AVX main loop over 32-float chunks implicitly (loop of 8)
// and supports the masked tail (partial) via `count`. When count < pushed, the masked
// lanes of the group are handled by the caller with bait; here we simply clamp.
static inline void rotor_batch(const float* src, float* sin_out, float* cos_out, size_t n) {
    for (size_t i = 0; i < n; ++i) {
        float si, co;
        rotor_sincos(src[i], si, co);
        sin_out[i] = si;
        cos_out[i] = co;
    }
}

} // namespace detkernel