#include "vlog.hpp"
#include <cmath>
#include <cstring>

namespace vlog {

namespace {

// Float bit patterns used by the AVX kernel's integer lane ops.
constexpr uint32_t C_2_OVER_3_BITS = 0x3f2aaaab; // bits of 2/3f (0.6666666865f)
constexpr uint32_t C_MANTISSA_MASK = 0x007fffff; // 2^23 - 1 (low 23 mantissa bits)

// ln(2) as a single float constant (0x3f317218).
constexpr float C_LN2 = 0.6931471824645996f;

// Minimax polynomial coefficients for ln(1+x) on the reduced interval
// x in [-1/3, 1/3), read straight from the binary:
//   ln(1+x) ~= x + x^2 * (c7 + c6*x + c5*x^2 + c4*x^3 + c3*x^4 + c2*x^5 + c1*x^6)
constexpr float C_P1 = -0.15177205204963684f; // 0x181f82040
constexpr float C_P2 =  0.16964881122112274f; // 0x181f82020
constexpr float C_P3 = -0.16462457180023193f; // 0x181f82000
constexpr float C_P4 =  0.19822503626346588f; // 0x181f81fe0
constexpr float C_P5 = -0.25004664063453674f; // 0x181f81fc0
constexpr float C_P6 =  0.33336564898490906f; // 0x181f81fa0
constexpr float C_P7 = -0.5f;                 // 0x181f81f80

inline float fma_f(float a, float b, float c) { return __builtin_fmaf(a, b, c); }

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

// Fast path: input is a positive normal float -> ln(x).
//
// Range reduction via the "2/3" magic (equivalent to the kernel's
// vpsubd/vpsrad 0x17/vpand/vpaddd): the mantissa is folded into
// m in [2/3, 4/3) and the integer exponent e is recovered, so that
//   y = m * 2^e   =>   ln(y) = e*ln2 + ln(m),
// with ln(m) = ln(1+x), x = m-1 in [-1/3, 1/3), from the minimax polynomial.
inline float ln_fast(float y) {
    uint32_t b  = bit_u32(y);
    uint32_t t  = b - C_2_OVER_3_BITS;                      // vpsubd (wrapping)
    int32_t  e  = static_cast<int32_t>(t) >> 23;            // vpsrad 0x17 (exponent)
    uint32_t mb = (t & C_MANTISSA_MASK) + C_2_OVER_3_BITS;  // vpand + vpaddd
    float m = bit_f32(mb);
    float x = m - 1.0f;                                     // vsubps  (x in [-1/3, 1/3))
    float e_ = static_cast<float>(e);                       // vcvtdq2ps

    // Horner evaluation (mirrors the vfmadd231ps/vfmadd213ps chain).
    float p = C_P2;
    p = fma_f(x, C_P1, p);
    p = fma_f(p, x, C_P3);
    p = fma_f(p, x, C_P4);
    p = fma_f(p, x, C_P5);
    p = fma_f(p, x, C_P6);
    p = fma_f(p, x, C_P7);

    float q = x * p;          // x*P(x)
    q = fma_f(q, x, x);       // x^2*P(x) + x  ~= ln(1+x) = ln(m)
    return fma_f(e_, C_LN2, q); // e*ln2 + ln(m) == ln(y)
}

} // namespace

void log_f32(const float* src, float* dst, uint32_t n) {
    for (uint32_t i = 0; i < n; ++i) {
        float y = src[i];
        uint32_t b = bit_u32(y);

        // Fast path iff the lane is a positive normal float:
        // (int32)(b + 0x00800000) >= 0x01000000  <=>  b in [0x00800000, 0x7f7fffff].
        if (b >= 0x00800000u && b <= 0x7f7fffffu) {
            dst[i] = ln_fast(y);
        } else {
            // Slow path: zero/denormal/negative/Inf/NaN. The kernel dispatches to
            // a scalar double-precision Cody-Waite ln (0x1802a2fc0); std::log is the
            // structurally equivalent reference for these edge inputs.
            dst[i] = static_cast<float>(std::log(static_cast<double>(y)));
        }
    }
}

} // namespace vlog