g_math 0.4.25

Multi-domain fixed-point arithmetic with geometric extension: Lie groups, manifolds, ODE solvers, tensors, fiber bundles — pure Rust, zero-float, deterministic
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137

#!/usr/bin/env python3
"""Generate gMath decimal contract-validation golden tables.

Independent mpmath oracle for gMath's DecimalFixed<SCALE> domain. Emits the
same mode-agnostic format used by mootable/decimal-scaled so a single Rust
grader handles both this (oracle #1) and that corpus (oracle #2):

    <input_raw>\t<floor_raw>\t<cls>

  input_raw  storage integer; mathematical value x = input_raw / 10**SCALE.
  floor_raw  floor(f(x) * 10**SCALE), toward negative infinity (sign here).
  cls        fractional class of f(x)*10**SCALE - floor_raw, in [0, 1):
               Z exact, L below-half, E exact-half-tie, G above-half.

NO floats cross the boundary: mpmath computes at high guard precision and we
emit exact integer strings. The Rust side is pure integer arithmetic.

Deterministic: fixed seed, sorted inputs -> byte-identical reruns.

Usage:  python3 scripts/gen_decimal_golden.py
"""
from __future__ import annotations

import random
from pathlib import Path

from mpmath import mp, mpf, floor as mp_floor, exp, ln, sqrt, sin

# Scale 28 lands on decimal-scaled's d38_s28 grid point (both i128), so the
# phase-2 second-oracle comparison is apples-to-apples at the same grid cell.
SCALE = 28
SCALE_POW = 10 ** SCALE
I128_MAX = (1 << 127) - 1
# Stay clear of the i128 edge so neither input_raw nor floor_raw overflows.
RANGE_LIMIT = mpf(I128_MAX) / mpf(SCALE_POW) * mpf("0.95")  # ~1.6e10

# Guard precision: far wider than the storage scale so L/E/G is unambiguous.
mp.dps = 2 * SCALE + 80

OUT_DIR = Path(__file__).resolve().parent.parent / "tests" / "oracle_golden"
PER_FUNC = 64


def floor_and_class(value: mpf) -> tuple[int, str] | None:
    """(floor_raw, cls) for value at SCALE, or None if it leaves i128 range."""
    scaled = value * mpf(SCALE_POW)
    if abs(scaled) >= mpf(I128_MAX):
        return None
    floor_int = int(mp_floor(scaled))
    frac = scaled - mpf(floor_int)
    half = mpf("0.5")
    if frac == 0:
        cls = "Z"
    elif frac < half:
        cls = "L"
    elif frac == half:
        cls = "E"
    else:
        cls = "G"
    return floor_int, cls


def to_raw(x: mpf) -> int:
    """Storage integer for input value x (round to nearest, ties away)."""
    scaled = x * mpf(SCALE_POW)
    floor_int = int(mp_floor(scaled))
    frac = scaled - mpf(floor_int)
    return floor_int + 1 if frac >= mpf("0.5") else floor_int


def gen_inputs(name: str, rng: random.Random) -> list[mpf]:
    """In-range, well-conditioned input values for one function."""
    xs: list[mpf] = []
    if name == "exp":
        # output exp(x) must stay < RANGE_LIMIT  ->  x < ln(RANGE_LIMIT) ~ 23.5
        hi = float(ln(RANGE_LIMIT))
        xs += [mpf(rng.uniform(-hi, hi)) for _ in range(PER_FUNC - 4)]
        xs += [mpf(0), mpf("0.5"), mpf(-1), mpf(1)]
    elif name == "ln":
        # x in (0, RANGE_LIMIT); output ln(x) within range for x > ~5e-9
        lo, hi = -20.0, float(ln(RANGE_LIMIT))
        xs += [exp(mpf(rng.uniform(lo, hi))) for _ in range(PER_FUNC - 3)]
        xs += [mpf(1), mpf(2), mpf("0.5")]
    elif name == "sqrt":
        hi = float(RANGE_LIMIT)
        xs += [mpf(rng.uniform(0, hi)) for _ in range(PER_FUNC - 4)]
        # perfect squares -> exact (Z) coverage
        xs += [mpf(4), mpf(9), mpf("2.25"), mpf(0)]
    elif name == "sin":
        xs += [mpf(rng.uniform(-20.0, 20.0)) for _ in range(PER_FUNC - 3)]
        xs += [mpf(0), mpf(1), mpf(-1)]
    else:
        raise ValueError(name)
    return xs


FUNCS = {"exp": exp, "ln": ln, "sqrt": sqrt, "sin": sin}


def main() -> None:
    OUT_DIR.mkdir(parents=True, exist_ok=True)
    for name, fn in FUNCS.items():
        rng = random.Random(f"gmath-decimal-golden-{name}-s{SCALE}")
        rows: list[tuple[int, int, str]] = []
        seen: set[int] = set()
        for x in gen_inputs(name, rng):
            raw = to_raw(x)
            if raw in seen:
                continue
            seen.add(raw)
            xv = mpf(raw) / mpf(SCALE_POW)
            if name in ("ln", "sqrt") and xv < 0:
                continue
            if name == "ln" and xv == 0:
                continue
            res = floor_and_class(fn(xv))
            if res is None:
                continue  # out of representable range; skip
            floor_raw, cls = res
            rows.append((raw, floor_raw, cls))
        rows.sort()
        path = OUT_DIR / f"{name}_s{SCALE}.txt"
        with path.open("w") as f:
            f.write("# gMath decimal contract golden table (oracle #1)\n")
            f.write("# generated by scripts/gen_decimal_golden.py\n")
            f.write(f"# function={name}  scale={SCALE}  guard_dps={mp.dps}\n")
            f.write("# each line: <input_raw>\\t<floor_raw>\\t<cls>\n")
            f.write("# value x = input_raw / 10**scale; "
                    "floor_raw = floor(f(x)*10**scale) toward -inf;\n")
            f.write("# cls: Z exact, L below-half, E half-tie, G above-half.\n")
            for raw, floor_raw, cls in rows:
                f.write(f"{raw}\t{floor_raw}\t{cls}\n")
        print(f"{path.name}: {len(rows)} rows")


if __name__ == "__main__":
    main()