1import json, math, random
 2import numpy as np
 3
 4EPS = 1e-30
 5
 6def norm_pair(p):
 7    p = np.asarray(p, dtype=np.float64)
 8    s = max(float(np.max(np.abs(p))), EPS)
 9    return p / s
10
11def tensor(a, b):
12    return np.array([a[0]*b[0], a[1]*b[1]], dtype=np.float64)
13
14def direct_sum(a, b):
15    return np.array([a[0]*b[0], a[1]*b[0] + a[0]*b[1]], dtype=np.float64)
16
17def val(a):
18    if a[0] == 0: return math.inf
19    return a[1] / a[0]
20
21def enc_unsigned(x):
22    return norm_pair([1.0, max(float(x), 0.0)])
23
24def enc_signed(x):
25    return enc_unsigned(max(x, 0)), enc_unsigned(max(-x, 0))
26
27def signed_mul(x, y, normalize=True):
28    xp, xn = x; yp, yn = y
29    pp, nn = tensor(xp, yp), tensor(xn, yn)
30    pn, npair = tensor(xp, yn), tensor(xn, yp)
31    pos = direct_sum(pp, nn)
32    neg = direct_sum(pn, npair)
33    if normalize: pos, neg = norm_pair(pos), norm_pair(neg)
34    return pos, neg
35
36def signed_value(x):
37    return val(x[0]) - val(x[1])
38
39def algebra_check(seed=7, n=1000):
40    rng = np.random.default_rng(seed)
41    max_mul = max_add = max_gauge = 0.0
42    for _ in range(n):
43        a, b = 10**rng.uniform(-8, 8, 2)
44        A, B = enc_unsigned(a), enc_unsigned(b)
45        max_mul = max(max_mul, abs(val(norm_pair(tensor(A,B))) - a*b) / max(1, abs(a*b)))
46        max_add = max(max_add, abs(val(norm_pair(direct_sum(A,B))) - (a+b)) / max(1, abs(a+b)))
47        c = 10**rng.uniform(-10, 10)
48        max_gauge = max(max_gauge, abs(val(A*c)-a) / max(1, abs(a)))
49    max_signed = 0.0
50    for _ in range(n):
51        x, y = rng.uniform(-10, 10, 2)
52        got = signed_value(signed_mul(enc_signed(x), enc_signed(y)))
53        max_signed = max(max_signed, abs(got-x*y)/max(1,abs(x*y)))
54    return dict(max_unsigned_mul_relerr=max_mul, max_unsigned_add_relerr=max_add,
55                max_signed_mul_relerr=max_signed, max_gauge_relerr=max_gauge)
56
57def product_stress(dtype=np.float32, depth=200, factor=1.25):
58    # Positive products: normalized pair rails are bounded, but multiplication itself
59    # still performs the same potentially overflowing products as the scalar control.
60    f = np.array(factor, dtype=dtype)
61    scalar = np.array(1.0, dtype=dtype)
62    pair = np.array([1.0, factor], dtype=dtype)
63    finite_scalar = True; finite_pair = True; max_rail = 0.0
64    for _ in range(depth):
65        scalar = scalar * f
66        pair = pair * np.array([1.0, factor], dtype=dtype)
67        s = np.max(np.abs(pair))
68        if np.isfinite(s) and s > 0: pair = pair / s
69        finite_scalar &= bool(np.isfinite(scalar))
70        finite_pair &= bool(np.all(np.isfinite(pair)))
71        max_rail = max(max_rail, float(np.max(np.abs(pair))) if np.all(np.isfinite(pair)) else math.inf)
72    return dict(dtype=str(dtype), depth=depth, factor=factor,
73                scalar_finite=finite_scalar, pair_finite=finite_pair,
74                scalar_value=float(scalar) if np.isfinite(scalar) else None,
75                pair_rails=pair.tolist() if np.all(np.isfinite(pair)) else None,
76                decoded_pair=float(pair[1]/pair[0]) if pair[0] != 0 and np.all(np.isfinite(pair)) else None,
77                max_normalized_rail=max_rail)
78
79def cancellation_check():
80    # A signed representation with both sign rails active has cancellation; compare
81    # decoded error after repeated multiplication against ordinary scalar arithmetic.
82    vals = [1.0 + 1e-5, 1.0, 1.0 - 1e-5]
83    x = enc_signed(vals[0]); y = enc_signed(vals[1]); z = enc_signed(vals[2])
84    pair = signed_mul(signed_mul(x,y),z)
85    return dict(expected=math.prod(vals), decoded=signed_value(pair),
86                abs_error=abs(signed_value(pair)-math.prod(vals)))
87
88def main():
89    np.seterr(over='ignore', under='ignore', invalid='ignore', divide='ignore')
90    out = dict(algebra=algebra_check(), cancellation=cancellation_check(),
91               stress=[product_stress(np.float32, 200, 1.25),
92                       product_stress(np.float32, 100, 1.1),
93                       product_stress(np.float64, 1000, 1.1)])
94    print(json.dumps(out, indent=2, allow_nan=False))
95
96if __name__ == '__main__': main()