g_math 0.6.8

#!/usr/bin/env python3
"""Generate mpmath reference values for tests/data/compute_tier_refs.rs.

For each (profile, function, input): emit
  (x_compute_raw, expected_compute_raw, expected_storage_raw)
where compute raw = round(v * 2^CFB), storage raw = round(v * 2^FB).
All integers emitted as i128-compatible literals; inputs are chosen exactly
representable at Q16.16 (multiples of 2^-6), so both profiles share them.

The realtime profile gets one table per GMATH_FRAC_BITS (2..=30) and a `refs`
facade that selects the build's split. A split keeps an input only if the
input is exact at that split and both the input and the result fit its
storage range. The FRAC_BITS=16 tables are the original q16_16 tables.
"""
from mpmath import mp, mpf, exp, ln, sqrt, sinh, cosh
from fractions import Fraction

mp.dps = 60


INPUTS = {
    "exp":      ["-3.5", "-0.75", "0.0625", "0.5", "1.25", "2", "3.5"],
    "ln":       ["0.0625", "0.5", "1.25", "2", "7.75"],
    "sqrt":     ["0.0625", "0.5", "2", "7.75", "100"],
    "sigmoid":  ["-7", "-2.5", "-0.5", "0", "0.5", "2.5", "7"],
    "softplus": ["-7", "-2.5", "-0.5", "0", "0.5", "2.5", "7"],
    "ln1p":     ["-0.9375", "-0.5", "0.0625", "1", "7.75"],
    "sinh":     ["-2.5", "-0.5", "0.5", "2.5"],
    "cosh":     ["-2.5", "-0.5", "0.5", "2.5"],
}

FUNCS = {
    "exp": exp,
    "ln": ln,
    "sqrt": sqrt,
    "sigmoid": lambda x: 1 / (1 + exp(-x)),
    "softplus": lambda x: ln(1 + exp(x)),
    "ln1p": lambda x: ln(1 + x),
    "sinh": sinh,
    "cosh": cosh,
}

def to_raw(v, bits):
    """round-to-nearest integer of v * 2^bits (ties away from zero)."""
    scaled = v * (mpf(2) ** bits)
    from mpmath import floor
    f = floor(scaled)
    frac = scaled - f
    n = int(f)
    if frac > mpf("0.5"):
        n += 1
    elif frac == mpf("0.5"):
        n += 1 if n >= 0 else 0  # ties toward +inf, matching round-bit add
    return n

def rows_for(name, fb, cfb, storage_bits):
    limit = mpf(2) ** (storage_bits - 1 - fb)
    rows = []
    for xs_s in INPUTS[name]:
        xfrac = Fraction(xs_s)
        # input raw at storage then exact at compute (input must be exactly representable)
        if (xfrac * (2 ** fb)).denominator != 1:
            continue
        v = FUNCS[name](mpf(xs_s))
        if abs(mpf(xs_s)) >= limit or abs(v) >= limit - 1:
            continue
        x_compute_raw = int(xfrac * (2 ** cfb))
        rows.append((x_compute_raw, to_raw(v, cfb), to_raw(v, fb)))
    return rows

def emit_tables(out, fb, cfb, storage_bits, indent):
    for name in INPUTS:
        out.append(f"{indent}pub const {name.upper()}: &[(i128, i128, i64)] = &[")
        for r in rows_for(name, fb, cfb, storage_bits):
            out.append(f"{indent}    ({r[0]}, {r[1]}, {r[2]}),")
        out.append(f"{indent}];")

out = []
out.append("// GENERATED by scripts/generate_compute_tier_refs.py — mpmath 60-digit references.")
out.append("// (x_compute_raw, expected_compute_raw, expected_storage_raw)")
out.append("// compute raw = round(v * 2^COMPUTE_FRAC_BITS); storage raw = round(v * 2^FRAC_BITS).")
out.append("")

# realtime: one module per GMATH_FRAC_BITS, and a facade for the build's split
SPLITS = list(range(2, 31))
for fb in SPLITS:
    out.append('#[cfg(table_format = "q16_16")]')
    out.append(f"pub mod refs_f{fb} {{")
    emit_tables(out, fb, 2 * fb, 32, "    ")
    out.append("}")
    out.append("")
out.append('#[cfg(table_format = "q16_16")]')
out.append("pub mod refs {")
for name in INPUTS:
    arms = ", ".join(f"{fb} => super::refs_f{fb}::{name.upper()}" for fb in SPLITS if fb != 16)
    out.append(f"    pub const {name.upper()}: &[(i128, i128, i64)] =")
    out.append(f"        match g_math::fixed_point::frac_config::FRAC_BITS {{ {arms}, _ => super::refs_f16::{name.upper()} }};")
out.append("}")
out.append("")

out.append('#[cfg(table_format = "q32_32")]')
out.append("pub mod refs {")
emit_tables(out, 32, 64, 64, "    ")
out.append("}")
out.append("")

print("\n".join(out))