g_math 0.4.32

Multi-domain fixed-point arithmetic with geometric extension: Lie groups, manifolds, ODE solvers, tensors, fiber bundles — pure Rust, zero-float, deterministic
Documentation
#!/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.
"""
from mpmath import mp, mpf, exp, ln, sqrt, sinh, cosh
from fractions import Fraction

mp.dps = 60

PROFILES = {"q16_16": (16, 32), "q32_32": (32, 64)}

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

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("")
for prof, (fb, cfb) in PROFILES.items():
    out.append(f'#[cfg(table_format = "{prof}")]')
    out.append("pub mod refs {")
    for name, xs in INPUTS.items():
        rows = []
        for xs_s in xs:
            xfrac = Fraction(xs_s)
            # input raw at storage then exact at compute (input is exactly representable)
            x_storage_raw = xfrac * (2 ** fb)
            assert x_storage_raw.denominator == 1, f"input {xs_s} not exact at Q{fb}"
            x_compute_raw = int(xfrac * (2 ** cfb))
            v = FUNCS[name](mpf(xs_s))
            rows.append((x_compute_raw, to_raw(v, cfb), to_raw(v, fb)))
        out.append(f"    pub const {name.upper()}: &[(i128, i128, i64)] = &[")
        for r in rows:
            out.append(f"        ({r[0]}, {r[1]}, {r[2]}),")
        out.append("    ];")
    out.append("}")
    out.append("")

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