# * Soft Math Library for Viper
# * Implements a self-contained software floating-point math runtime for float (binary32) and double (binary64),
# * without external libm dependencies. Designed for embedded, bare-metal, or custom VM environments where hardware
# * FPU math routines or system-provided standard math functions are unavailable, undesirable, or require fully
# * auditable source.
# *
# * Scope and feature set:
# * - Fundamental mathematical constants: π, e, π/2, π/4, reciprocal π, 2/π, ln(2), ln(10), 2/sqrt(π), DBL_MAX.
# * - Angle conversion: radians ↔ degrees helpers.
# * - Trigonometric functions: sin, cos, tan, sinf, cosf, tanf. Argument reduction is O(1): a three-constant
# *   Cody-Waite step for |x| < 2^20, and a Payne-Hanek step (20-word table of 2/π) above. Both yield r ∈ [-π/4, π/4]
# *   as a double-double fed to the truncated Taylor series.
# * - Inverse trigonometric functions: asin, acos, atan, atan2. asin uses a minimax polynomial suitable for [-1, 1];
# *   atan relies on a Taylor core limited to |x| ≤ 1 combined with range folding via reciprocal and quadrant logic.
# * - Hyperbolic trigonometry: sinh, cosh, tanh, built atop exp and exp(-x).
# * - Exponential and logarithm: exp uses range reduction into [-ln2/2, ln2/2] followed by a Taylor polynomial for the
# *   reduced residual, then scales by powers of two. Natural log implements a binary normalization pass to bring the
# *   mantissa into [1, 2), then uses the rational substitution y=(s-1)/(s+1) with a convergent series.
# * - Derived logarithms and exponentials: log10, log2, exp2, expm1, log1p.
# * - General power: pow for real bases and exponents, implemented as exp(y * log(x)), with special-case handling for
# *   zero and negative bases. cbrt is a wrapper over pow using the 1/3 exponent.
# * - Root and distance: Newton–Raphson sqrt / sqrtf with fixed 10 iterations; hypot for Euclidean norm.
# * - Rounding primitives: floor, ceil, round, trunc, with float variants floorf, ceilf.
# * - Modulo: fmod, fmodf following the C remainder convention where the result carries the sign of the dividend.
# * - Decomposition: modf to split a double into integer and fractional parts via a pointer output.
# * - Inverse hyperbolics: asinh, acosh, atanh expressed through log and sqrt compositions.
# * - Special functions: Gamma using Stirling’s approximation; error function erf / complementary erfc using a simple
# *   rational exponential approximation.
# * - Combinatorics: factorial, binomial combination, permutation, returning double for wider dynamic range.
# * - Integer helpers: abs for int, labs for long, fabs / fabsf for floating-point absolute value.
# * - IEEE-754 predicate helpers: isnan and isinf implemented via type punning union to inspect binary64 bit patterns.
# *   These avoid numeric conversion which would destroy NaN/inf payloads; they examine exponent and mantissa bits directly.
# *
# * Precision notes:
# * All polynomial cores are truncated series, not infinite expansions. Errors are bounded for the reduced interval but
# * will accumulate through composition (e.g. pow depends on log and exp; trigonometric inverses depend on asin/atan cores).
# * This implementation prioritizes determinism and source transparency over last-bit IEEE compliance. It is not a drop-in
# * replacement for a high-quality libm with accurate argument reduction, correct rounding, and exhaustive edge-case
# * handling. Input values outside the designed convergence domains will see degraded accuracy.
# *
# * Range and edge behavior:
# * exp clamps extreme inputs: values above approximately 709 return DBL_MAX, values below approximately -745 return 0.0,
# * preventing uncontrolled overflow in downstream code such as softmax. log returns NEG_1E308 for non-positive inputs
# * as a sentinel instead of NaN. Division-by-zero conditions in tan/tanf return zero as a simple fallback; other edge
# * cases propagate NaN where the underlying computation produces it. The isnan/isinf predicates use direct bit inspection
# * and respect the binary64 encoding rules for quiet/signaling NaN and signed infinities.
# *
# * Usage notes:
# * This library assumes Viper’s double maps to IEEE 754 binary64 and float maps to IEEE 754 binary32. Type punning in
# * isnan/isinf relies on a union reinterpretation of the underlying memory layout; this is portable only on targets
# * with standard little-endian or big-endian IEEE float storage and must be reviewed for exotic memory models.
# * Polynomial loops are manually folded into fixed for-loops over coefficient arrays to remove redundant source lines
# * while preserving the original floating-point evaluation order, which matters for reproducibility and error accumulation.
# *
# * License: MIT


define U_M_PI:   double = 3.14159265358979323846
define U_M_E:    double = 2.71828182845904523536
define U_M_PI_2: double = 1.57079632679489661923
define U_M_PI_4: double = 0.78539816339744830962
define U_M_1_PI: double = 0.31830988618379067154
define U_M_2_PI: double = 0.63661977236758134308

define U_M_LN2:  double = 0.69314718055994530942
define U_M_LN10: double = 2.30258509299404568402

define U_M_2_SQRT_PI: double = 1.77245385090551602730

define DBL_MAX: double = 1.7976931348623157e308
define NEG_1E308: double = -1e308

# Argument reduction constants for sin/cos/tan: pi/2 split into HI+MID+LO and 2/pi (all exact doubles),
# plus the 2^20 limit below which Cody-Waite applies (so each fn*HK product stays exact in 53 bits).
define U_M_PIO2_HI:  double = 1.57079632673412561417e+00
define U_M_PIO2_MID: double = 6.07710050650619224932e-11
define U_M_PIO2_LO:  double = 3.52155986518320018165e-27
define U_M_INV_PIO2: double = 6.36619772367581382433e-01
define U_M_REDUCE_LIMIT: double = 1048576.0

def radians(degrees: double) -> double: return degrees * (U_M_PI / 180.0)
def degrees(radians: double) -> double: return radians * (180.0 / U_M_PI)

def _sin_core(x: double) -> double:
    """sin(r) Taylor series, 8 terms, valid on |r| <= pi/4."""
    x2: double = x * x
    result: double = x
    term: double = x
    sign: double = -1.0
    coeffs: const [8]double = [6.0, 20.0, 42.0, 72.0, 110.0, 156.0, 210.0, 272.0]
    for c in coeffs:
        term = term * x2 / c
        result += sign * term
        sign = -sign
    return result


def _cos_core(x: double) -> double:
    """cos(r) Taylor series, 8 terms, valid on |r| <= pi/4."""
    x2: double = x * x
    result: double = 1.0
    term:   double = 1.0
    sign:   double = -1.0
    coeffs: const [8]double = [2.0, 12.0, 30.0, 56.0, 90.0, 132.0, 182.0, 240.0]
    for c in coeffs:
        term = term * x2 / c
        result += sign * term
        sign = -sign
    return result


def _mul_add(a: uint64, b: uint64, c: uint64, lo: uint64*, hi: uint64*):
    """Write (hi:lo) = a*b + c using 32-bit limbs for the full 128-bit product."""
    a0: uint64 = a & 0xFFFFFFFF
    a1: uint64 = a >> 32
    b0: uint64 = b & 0xFFFFFFFF
    b1: uint64 = b >> 32
    w0: uint64 = a0 * b0
    t: uint64 = a1 * b0 + (w0 >> 32)
    w1: uint64 = t & 0xFFFFFFFF
    w2: uint64 = t >> 32
    w1 = a0 * b1 + w1
    h: uint64 = a1 * b1 + w2 + (w1 >> 32)
    l: uint64 = (w1 << 32) | (w0 & 0xFFFFFFFF)
    l2: uint64 = l + c
    if l2 < l: h = h + 1
    *lo = l2
    *hi = h


def _two_prod(a: double, b: double, ph: double*, pl: double*):
    """Write a*b as the exact double-double (ph + pl) via Dekker splitting (no FMA required)."""
    split: double = 134217729.0
    p: double = a * b
    ca: double = split * a
    ahi: double = ca - (ca - a)
    alo: double = a - ahi
    cb: double = split * b
    bhi: double = cb - (cb - b)
    blo: double = b - bhi
    *ph = p
    *pl = ((ahi * bhi - p) + ahi * blo + alo * bhi) + alo * blo


def _rem_pio2(x: double, rh: double*, rl: double*) -> int:
    """O(1) reduction of finite x >= 0: x = n*(pi/2) + r with |r| <= pi/4; returns n mod 4, writes r = rh + rl."""
    if x < U_M_REDUCE_LIMIT:
        # Cody-Waite three-step with pi/2 = HI + MID + LO; each fn*HK product is exact for fn < 2^20
        q: double = x * U_M_INV_PIO2
        k: int = int(q)
        if q - double(k) > 0.5: k = k + 1
        fn: double = double(k)
        r: double = x - fn * U_M_PIO2_HI
        w1: double = fn * U_M_PIO2_MID
        y0: double = r - w1
        w2: double = fn * U_M_PIO2_LO
        rhv: double = y0 - w2
        *rh = rhv
        *rl = ((r - y0) - w1) + ((y0 - rhv) - w2)
        return k & 3
    # Payne-Hanek: a 20-word (1280-bit) 2/pi table; multiply the 53-bit mantissa and read off the fraction
    twopi: const [20]uint64 = [
        0xA2F9836E4E441529, 0xFC2757D1F534DDC0, 0xDB6295993C439041, 0xFE5163ABDEBBC561,
        0xB7246E3A424DD2E0, 0x06492EEA09D1921C, 0xFE1DEB1CB129A73E, 0xE88235F52EBB4484,
        0xE99C7026B45F7E41, 0x3991D639835339F4, 0x9C845F8BBDF9283B, 0x1FF897FFDE05980F,
        0xEF2F118B5A0A6D1F, 0x6D367ECF27CB09B7, 0x4F463F669E5FEA2D, 0x7527BAC7EBE5F17B,
        0x3D0739F78A5292EA, 0x6BFB5FB11F8D5D08, 0x56033046FC7B6BAB, 0xF0CFBC209AF4361D
    ]
    xv: double = x
    xb: uint64 = 0
    *&xb = *&xv
    ef: uint64 = (xb >> 52) & 0x7FF
    m: uint64 = (xb & 0xFFFFFFFFFFFFF) | 0x10000000000000
    e: int = int(ef) - 1075
    shift: int = 1280 - e
    qw: int = shift >> 6
    s: int = shift & 63
    b: [21]uint64
    carry: uint64 = 0
    for j in range(20):
        lo: uint64 = 0
        hi: uint64 = 0
        _mul_add(m, twopi[19 - j], carry, &lo, &hi)
        b[j] = lo
        carry = hi
    b[20] = carry
    lo2: uint64 = b[qw] >> s
    if s == 63: lo2 = (b[qw] >> 63) | (b[qw + 1] << 1)
    k: int = int(lo2 & 3)
    a1: uint64 = b[qw] & ((uint64(1) << s) - 1)
    a2: uint64 = b[qw - 1]
    a3: uint64 = b[qw - 2]
    fhi: uint64 = a2
    flo: uint64 = a3
    if s != 0:
        sh: int = 64 - s
        fhi = (a1 << sh) | (a2 >> s)
        flo = (a2 << sh) | (a3 >> s)
    # fractional part f = fhi*2^-64 + flo*2^-128; split off the top 53 bits so df is exact
    f53: uint64 = fhi >> 11
    df: double = double(f53)
    de: double = double(fhi - (f53 << 11))
    fh: double = df * 1.1102230246251565e-16
    fl: double = de * 5.421010862427522e-20 + double(flo) * 2.9387358770557188e-39
    r0h: double = fh
    r0l: double = fl
    n: int = k
    if fh > 0.5 or (fh == 0.5 and fl > 0.0):
        n = (k + 1) & 3
        r0h = fh - 1.0
    ph: double = 0.0
    pl: double = 0.0
    _two_prod(r0h, U_M_PIO2_HI, &ph, &pl)
    *rh = ph
    *rl = pl + r0h * U_M_PIO2_MID + r0h * U_M_PIO2_LO + r0l * U_M_PIO2_HI
    return n


def sin(x: double) -> double:
    u: uint64 = 0
    xv: double = x
    *&u = *&xv
    ef: uint64 = (u >> 52) & 0x7FF
    if ef == 0x7FF: return x - x    # inf/nan -> NaN
    neg: int = int(u >> 63)
    ax: double = x if neg == 0 else -x
    rh: double = 0.0
    rl: double = 0.0
    n: int = _rem_pio2(ax, &rh, &rl)
    sr: double = _sin_core(rh) + _cos_core(rh) * rl
    cr: double = _cos_core(rh) - _sin_core(rh) * rl
    res: double = 0.0
    if   n == 0: res = sr
    elif n == 1: res = cr
    elif n == 2: res = -sr
    else:        res = -cr
    return -res if neg != 0 else res


def cos(x: double) -> double:
    u: uint64 = 0
    xv: double = x
    *&u = *&xv
    ef: uint64 = (u >> 52) & 0x7FF
    if ef == 0x7FF: return x - x    # inf/nan -> NaN
    neg: int = int(u >> 63)
    ax: double = x if neg == 0 else -x
    rh: double = 0.0
    rl: double = 0.0
    n: int = _rem_pio2(ax, &rh, &rl)
    sr: double = _sin_core(rh) + _cos_core(rh) * rl
    cr: double = _cos_core(rh) - _sin_core(rh) * rl
    if   n == 0: return cr
    elif n == 1: return -sr
    elif n == 2: return -cr
    else:        return sr


def tan(x: double) -> double: return sin(x) / cos(x)


def asin(x: double) -> double:
    if x > 1.0:  x = 1.0
    if x < -1.0: x = -1.0
    x2: double = x * x
    result: double = x
    result += 1.0  * x * x2                     / 6.0
    result += 3.0  * x * x2 * x2                / 40.0
    result += 5.0  * x * x2 * x2 * x2           / 112.0
    result += 35.0 * x * x2 * x2 * x2 * x2      / 1152.0
    result += 63.0 * x * x2 * x2 * x2 * x2 * x2 / 2816.0
    return result


def acos(x: double) -> double: return U_M_PI / 2.0 - asin(x)


def _atan_core(x: double) -> double:
    x2: double = x * x
    result: double = x
    term: double = x
    sign: double = -1.0
    coeffs: const [10]double = [3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0, 19.0, 21.0]
    for c in coeffs:
        term = term * x2 / c
        result += sign * term
        sign = -sign
    return result


def atan(x: double) -> double:
    ax: double = x
    if ax < 0.0:
        ax = -ax
    if ax > 1.0:
        return 0.0 - _atan_core(1.0 / ax) + U_M_PI_2 if x > 0.0 else _atan_core(1.0 / ax) - U_M_PI_2
    return _atan_core(x) if x >= 0.0 else 0.0 - _atan_core(ax)


def atan2(y: double, x: double) -> double:
    if   x > 0.0: return atan(y / x)
    elif x < 0.0:
        if y >= 0.0: return atan(y / x) + U_M_PI
        else:        return atan(y / x) - U_M_PI
    else:
        if   y > 0.0: return U_M_PI_2
        elif y < 0.0: return 0.0 - U_M_PI_2
        else:         return 0.0


def sinh(x: double) -> double:
    ex:     double = exp(x)
    e_negx: double = exp(-x)
    return (ex - e_negx) / 2.0


def cosh(x: double) -> double:
    ex:     double = exp(x)
    e_negx: double = exp(-x)
    return (ex + e_negx) / 2.0


def tanh(x: double) -> double:
    ex:     double = exp(x)
    e_negx: double = exp(-x)
    return (ex - e_negx) / (ex + e_negx)


def exp(x: double) -> double:
    """e^x
    Range reduction   Short Taylor: Directly expanding x with 18 terms of Taylor series diverges when |x| is slightly large
    (x=-10 gives 54.5, true value 4.5e-5; x=-30 gives 3.7e10), on the negative side it gives huge positive values,
    which can blow up the denominator in softmax → all probabilities ≈ 0. Instead, first do
    x = k*ln2 + r (|r| <= ln2/2), then calculate exp(r) and multiply by 2^k.
    """
    # First do a "finite and within range" positive check: comparisons with NaN are always false and would skip
    # the following two checks, and int(NaN) on x86 becomes INT_MIN -> the later 2^k scaling loop would run
    # 2.1 billion times, appearing as "frozen". Therefore NaN/overflow must be blocked here.
    if x <= 709.0 and x >= -745.0:
        pass
    else:
        if x > 709.0:  return DBL_MAX
        if x < -745.0: return 0.0
        return x # NaN: Return as is, let NaN propagate instead of crashing
    # Reduction: Fold x into [-ln2/2, ln2/2], accumulate the integer part into k (inline rounding to avoid forward reference)
    ln2: double = U_M_LN2
    q:   double = x / ln2
    k:   int = int(q)
    frac: double = q - double(k)
    if frac > 0.5: k = k + 1
    elif frac < -0.5: k = k - 1
    r: double = x - double(k) * ln2
    # IE # exp(r), |r| <= 0.3466: after 12 terms, truncation error < 1e-17
    result: double = 1.0
    term:   double = 1.0
    n: int = 1
    while n <= 12:
        term = term * r / double(n)
        result = result + term
        n += 1
    while k > 0:
        result = result * 2.0
        k = k - 1
    while k < 0:
        result = result * 0.5
        k = k + 1
    return result


def log(x: double) -> double:
    if x <= 0.0: return NEG_1E308
    exponent: int = 0
    while x >= 2.0:
        x /= 2.0
        exponent += 1
    while x < 1.0:
        x *= 2.0
        exponent -= 1
    y: double = (x - 1.0) / (x + 1.0)
    y2: double = y * y
    result: double = 2.0 * y
    term: double = 2.0 * y
    coeffs: const [10]double = [3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0, 19.0, 21.0]
    for c in coeffs:
        term *= y2
        result += term / c
    result += double(exponent) * U_M_LN2
    return result


def sqrt(x: double) -> double:
    if x < 0.0:  return 0.0
    if x == 0.0: return 0.0
    guess: double = x / 2.0
    for _ in range(10):
        guess = (guess + x / guess) / 2.0
    return guess


def abs(x: int) -> int: return x if x >= 0 else -x
def fabs(x: double) -> double: return x if x >= 0.0 else -x
def labs(x: long) -> long: return x if x >= 0 else -x

def factorial(n: int) -> double:
    if n < 0: return 0.0
    if n == 0 or n == 1: return 1.0
    result: double = 1.0
    for i in range(2, int(n + 1)):
        result *= i
    return result


def combination(n: int, k: int) -> double:
    if k < 0 or k > n: return 0.0
    if k == 0 or k == n: return 1.0
    k = k if k < n - k else n - k
    result: double = 1.0
    for i in range(1, int(k + 1)):
        result = result * (n - k + i) / i
    return result


def permutation(n: int, k: int) -> double:
    if k < 0 or k > n: return 0.0
    result: double = 1.0
    for i in range(0, int(k)):
        result *= (n - i)
    return result

def pow(x: double, y: double) -> double:
    if x == 0.0: return 0.0 if y > 0.0 else 1.0
    if x  < 0.0: return exp(y * log(-x))
    return exp(y * log(x))


def  powf(x: double, y: double) -> double: return pow(x, y)
def  cbrt(x: double           ) -> double: return pow(x, 1.0 / 3.0)
def hypot(x: double, y: double) -> double: return sqrt(x * x + y * y)


def floor(x: double) -> double:
    if x >= 0.0: return double(int(x))
    i: int = int(x)
    if x < double(i): return double(i - 1)
    else:             return double(i)


def ceil(x: double) -> double:
    if x < 0.0: return double(int(x))
    i: int = int(x)
    if x > double(i): return double(i + 1)
    else:             return double(i)


def round(x: double) -> double: return floor(x + 0.5)
def trunc(x: double) -> double: return floor(x) if x >= 0.0 else ceil(x)


def fmod(x: double, y: double) -> double:
    if y == 0.0: return 0.0
    result: double = x - y * floor(x / y)
    if result < 0.0: result += y
    return result


def fmodf(x: float, y: float) -> float:
    if y == 0.0: return 0.0
    result: float = x - y * floorf(x / y)
    if result < 0.0: result += y
    return result


def modf(x: double, iptr: double*) -> double:
    int_part: double = floor(x)
    *iptr = int_part
    return x - int_part


def log10(x: double) -> double: return log(x) / U_M_LN10
def  log2(x: double) -> double: return log(x) / U_M_LN2


def  exp2(x: double) -> double: return pow(2.0, x)
def expm1(x: double) -> double: return exp(x) - 1.0

def log1p(x: double) -> double: return log(1.0 + x)

def asinh(x: double) -> double: return log(x + sqrt(x * x + 1.0))
def acosh(x: double) -> double: return log(x + sqrt(x * x - 1.0))
def atanh(x: double) -> double: return (log(1.0 + x) - log(1.0 - x)) / 2.0


def gamma(x: double) -> double:
    if x <= 0.0: return 0.0
    if x == 1.0 or x == 2.0: return 1.0
    ln_gamma: double = 0.5 * log(2.0 * U_M_PI / x) + x * log(x / U_M_E)
    return exp(ln_gamma)


def erf(x: double) -> double:
    sign: double = 1.0 if x >= 0.0 else -1.0
    abs_x: double = fabs(x)
    t0: double = 4.0 * abs_x * abs_x / U_M_PI
    exp_term: double = exp(-t0)
    return sign * sqrt(1.0 - exp_term)


def erfc(x: double) -> double:
    return 1.0 - erf(x)

def sqrtf(x: float) -> float:
    if x <  0.0: return 0.0
    if x == 0.0: return 0.0
    guess: float = x / 2.0
    for i in range(10):
        guess = (guess + x / guess) / 2.0
    return guess

def sinf(x: float) -> float: return float(sin(double(x)))
def cosf(x: float) -> float: return float(cos(double(x)))

def tanf(x: float) -> float:
    c: float = cosf(x)
    if c == 0.0: return 0.0
    return sinf(x) / c

def fabsf(x: float) -> float: return x if x >= 0.0 else -x

def floorf(x: float) -> float:
    if x >= 0.0: return float(int(x))
    i: float = float(int(x))
    if x < i: return i - 1
    else:     return i

def ceilf(x: float) -> float:
    if x < 0.0: return float(int(x))
    i: float = float(int(x))
    if x > i: return i + 1
    else: return i


union _U:
    d: double
    u: uint64

def isnan(x: double) -> int:
    # IEEE-754 bit reinterpretation (can't use uint64(x) numerical conversion: that's fptosi, which loses the bit pattern completely
    # → exp_part is always 0 → always returns 0). Bitwise copy
    xv: double = x
    u: uint64 = 0
    *&u = *&xv
    exp_part: uint64 = (u >> 52) & 0x7FF
    mantissa_part: uint64 = u & 0xFFFFFFFFFFFFF
    return 1 if exp_part == 0x7FF and mantissa_part != 0 else 0

def isinf(x: double) -> int:
    xv2: double = x
    u2: uint64 = 0
    *&u2 = *&xv2
    exp_part2: uint64 = (u2 >> 52) & 0x7FF
    mantissa_part2: uint64 = u2 & 0xFFFFFFFFFFFFF
    return 1 if exp_part2 == 0x7FF and mantissa_part2 == 0 else 0