g_math 0.6.8

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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182

#!/usr/bin/env python3
"""References for the 0.6.7 wide-tier gate functions (tests/data/gate_refs.rs).

Q64.64 functions (`g_math::wide`): sigmoid, softplus, silu, sqrt. Inputs are
Q.40 raws (the value k / 2^40), given to the function as `k << 24`. For each,
the reference is the exact value rounded to nearest (ties toward +infinity)
at 64, 40 and 20 fractional bits, from mpmath at 140 digits (sqrt: exact
integer arithmetic), so the test can measure the Q64.64 error in units and
check the narrowed results for exact equality.

Storage functions (`fused::sigmoid_mul`, `fused::entropy`): dyadic inputs
with at most 8 fraction bits, stored exactly at every gated split; the
reference is a 120-digit literal the test parses with the exact literal
parser, which rounds it to nearest at the build's split.

Integer functions (`wide::mul_div_*`, `wide::try_ratio_from_str`): exact
Python integers and `Fraction`.

Usage: python3 scripts/generate_gate_refs.py > tests/data/gate_refs.rs
"""
import random
from fractions import Fraction
from math import isqrt

from mpmath import mp, mpf, exp, log, floor

mp.dps = 140
SEED = 20261005
DIGITS = 120


def nearest(v, frac_bits):
    """round(v * 2^frac_bits), ties toward +infinity."""
    return int(floor(v * mpf(2) ** frac_bits + mpf(1) / 2))


def sigmoid(x):
    return 1 / (1 + exp(-x))


def softplus(x):
    return log(1 + exp(x))


def lit(v):
    return mp.nstr(v, DIGITS, strip_zeros=False, min_fixed=-10**9, max_fixed=10**9)


def dy_lit(fr):
    """Exact decimal literal of a dyadic Fraction."""
    sign = "-" if fr < 0 else ""
    fr = abs(fr)
    d = 0
    while (fr * 10 ** d).denominator != 1:
        d += 1
    n = int(fr * 10 ** d)
    if d == 0:
        return f"{sign}{n}"
    s = str(n).rjust(d + 1, "0")
    return f"{sign}{s[:-d]}.{s[-d:]}"


out = ["// Generated by scripts/generate_gate_refs.py. Do not edit.", ""]
rng = random.Random(SEED)

# ---- Q64.64 functions on Q.40 inputs -----------------------------------------
def q40_inputs(n, lo, hi):
    xs = [0, 1, -1, 1 << 40, -(1 << 40), 12 << 40, -(12 << 40), 39 << 40, -(39 << 40)]
    while len(xs) < n:
        mag = rng.choice([1, 4, 12, 40])
        xs.append(rng.randrange(int(lo * mag / 40 * 2 ** 40), int(hi * mag / 40 * 2 ** 40)))
    return [x for x in xs if lo * 2 ** 40 <= x <= hi * 2 ** 40]

for name, f in (("SIGMOID", sigmoid), ("SOFTPLUS", softplus), ("SILU", lambda x: x * sigmoid(x))):
    rows = []
    for k in q40_inputs(400, -40, 40):
        v = f(mpf(k) / mpf(2) ** 40)
        rows.append(f"({k}, {nearest(v, 64)}, {nearest(v, 40)}, {nearest(v, 20)})")
    out.append(f"/// (x at Q.40, f(x) at Q64.64, at Q.40, at Q.20), each the nearest value.")
    out.append(f"pub const {name}: &[(i128, i128, i128, i128)] = &[{', '.join(rows)}];")

# sqrt: exact. round(sqrt(x) * 2^f) for x = k / 2^40 is the nearest integer to sqrt(k * 2^(2f - 40)).
def sqrt_nearest(k, f):
    assert 2 * f >= 40
    n = k << (2 * f - 40)
    r = isqrt(n)
    return r + 1 if n > r * r + r else r

rows = []
ks = [0, 1, 2, 3, 1 << 40, 1 << 41, (1 << 40) + 1, 9 << 38, 1 << 86]
while len(ks) < 400:
    ks.append(rng.randrange(0, 1 << rng.choice([10, 30, 41, 50, 70, 86])))
for k in ks:
    rows.append(f"({k}, {sqrt_nearest(k, 64)}, {sqrt_nearest(k, 40)}, {sqrt_nearest(k, 20)})")
out.append("/// (x at Q.40, sqrt(x) at Q64.64, at Q.40, at Q.20), each the nearest value (exact integer arithmetic).")
out.append(f"pub const SQRT: &[(i128, i128, i128, i128)] = &[{', '.join(rows)}];")

# ---- fused::sigmoid_mul ---------------------------------------------------------
def dy(lo, hi, bits=8):
    return Fraction(rng.randrange(int(lo * 2 ** bits), int(hi * 2 ** bits)), 1 << bits)

def m_of(fr):
    return mpf(fr.numerator) / fr.denominator

rows = []
cases = [(Fraction(1), Fraction(0)), (Fraction(-1), Fraction(0)), (Fraction(3), Fraction(-100)), (Fraction(3), Fraction(100)),
         (Fraction(127), Fraction(1, 256)), (Fraction(-127), Fraction(-1, 256)), (Fraction(0), Fraction(5))]
while len(cases) < 240:
    cases.append((dy(-127, 127), dy(*rng.choice([(-2, 2), (-12, 12), (-60, 60)]))))
for x, g in cases:
    rows.append(f'("{dy_lit(x)}", "{dy_lit(g)}", "{lit(m_of(x) * sigmoid(m_of(g)))}")')
out.append("/// (x, gate, x * sigmoid(gate)); |x| < 128 so every gated split holds it.")
out.append(f"pub const SIGMOID_MUL: &[(&str, &str, &str)] = &[{', '.join(rows)}];")

# ---- fused::entropy ---------------------------------------------------------------
rows = []
for _ in range(60):
    n = rng.choice([1, 2, 5, 16, 64, 300])
    kind = rng.randrange(3)
    if kind == 0:      # a distribution on a 1/256 grid (sums to 1), zeros included
        cuts = sorted(rng.randrange(0, 257) for _ in range(n - 1))
        parts = [b - a for a, b in zip([0] + cuts, cuts + [256])]
        ws = [Fraction(p, 256) for p in parts]
    elif kind == 1:    # arbitrary small weights
        ws = [dy(0, 1) for _ in range(n)]
    else:              # weights above one (negative terms)
        ws = [dy(0, 4) for _ in range(min(n, 16))]
    h = -sum((m_of(w) * log(m_of(w)) for w in ws if w > 0), mpf(0))
    rows.append("(&[" + ", ".join(f'"{dy_lit(w)}"' for w in ws) + f'], "{lit(h)}")')
out.append("/// (weights, -sum(w ln w)).")
out.append(f"pub const ENTROPY: &[(&[&str], &str)] = &[{', '.join(rows)}];")

# ---- integer mul_div ---------------------------------------------------------------
I128_MIN, I128_MAX = -(1 << 127), (1 << 127) - 1
def r128():
    bits = rng.randrange(1, 128)
    v = rng.randrange(0, 1 << bits)
    return -v if rng.randrange(2) else v
rows = []
trip = [(I128_MAX, I128_MAX, I128_MAX), (I128_MIN, I128_MAX, I128_MIN), (I128_MIN, I128_MIN, I128_MAX), (7, -3, 2), (-7, 3, -2),
        (7, 3, 2), (-7, -3, -2), (1, 1, 2), (-1, 1, 2), (3, 1, 2), (-3, 1, 2), (0, 5, -9), (I128_MAX, 2, 1)]
while len(trip) < 600:
    trip.append((r128(), r128(), r128() or 1))
for a, b, d in trip:
    fl = (a * b) // d
    ne = (2 * a * b + d) // (2 * d)            # floor(p/d + 1/2): ties toward +infinity
    enc = lambda v: f"Some({v})" if I128_MIN <= v <= I128_MAX else "None"
    rows.append(f"({a}, {b}, {d}, {enc(fl)}, {enc(ne)})")
out.append("/// (a, b, d, floor(a*b/d), nearest(a*b/d) ties toward +infinity); None when it does not fit i128.")
out.append(f"pub const MUL_DIV: &[(i128, i128, i128, Option<i128>, Option<i128>)] = &[{', '.join(rows)}];")
rows = []
U = (1 << 128) - 1
ut = [(U, U, U), (U, U, 1), (U, 2, 3), (0, U, 7), (1 << 127, 2, 1), (1 << 127, 2, 2)]
while len(ut) < 300:
    ut.append((rng.randrange(0, 1 << rng.randrange(1, 129)), rng.randrange(0, 1 << rng.randrange(1, 129)), rng.randrange(1, 1 << rng.randrange(1, 129))))
for a, b, d in ut:
    q = (a * b) // d
    rows.append(f"({a}, {b}, {d}, {'Some(%d)' % q if q <= U else 'None'})")
out.append("/// (a, b, d, floor(a*b/d)); None when it does not fit u128.")
out.append(f"pub const MUL_DIV_U128: &[(u128, u128, u128, Option<u128>)] = &[{', '.join(rows)}];")

# ---- exact ratios -------------------------------------------------------------------
rows = []
lits = ["0.25", "10000000.0", "-1.5e1", "1e-5", ".5", "5.", "-0", "0.000", "+3", "1e3", "12.3400", "  7.5  ", "0.1", "1e38", "1e-38",
        "170141183460469231731687303715884105727", "0.5000000000000000000000000000000000000000000000000000", "1E2", "2.5e-1", "-0.125"]
for _ in range(180):
    s = rng.choice(["", "-", "+"]) + str(rng.randrange(0, 10 ** rng.randrange(1, 12)))
    if rng.randrange(3):
        s += "." + "".join(rng.choice("0123456789") for _ in range(rng.randrange(0, 14)))
    if rng.randrange(4) == 0:
        s += "e" + str(rng.randrange(-12, 13))
    lits.append(s)
for s in lits:
    body = s.strip().lstrip("+")
    fr = Fraction(body if body not in ("5.",) and not body.endswith(".") else body + "0") if not body.startswith(".") and not body.startswith("-.") else Fraction("0" + body)
    ok = I128_MIN <= fr.numerator <= I128_MAX and fr.denominator <= I128_MAX
    rows.append(f'("{s}", {"Some((%d, %d))" % (fr.numerator, fr.denominator) if ok else "None"})')
out.append("/// (literal, the exact fraction in lowest terms); None when it does not fit.")
out.append(f"pub const RATIO: &[(&str, Option<(i128, i128)>)] = &[{', '.join(rows)}];")
out.append('/// Strings outside the grammar.')
out.append('pub const RATIO_BAD: &[&str] = &["", " ", ".", "-", "e5", "1e", "1e+", "1.2.3", "0x10", "1/3", "pi", "1 2", "--1", "1e1.5", "٣"];')
print("\n".join(out))