import json, math, random
import numpy as np

EPS = 1e-30

def norm_pair(p):
    p = np.asarray(p, dtype=np.float64)
    s = max(float(np.max(np.abs(p))), EPS)
    return p / s

def tensor(a, b):
    return np.array([a[0]*b[0], a[1]*b[1]], dtype=np.float64)

def direct_sum(a, b):
    return np.array([a[0]*b[0], a[1]*b[0] + a[0]*b[1]], dtype=np.float64)

def val(a):
    if a[0] == 0: return math.inf
    return a[1] / a[0]

def enc_unsigned(x):
    return norm_pair([1.0, max(float(x), 0.0)])

def enc_signed(x):
    return enc_unsigned(max(x, 0)), enc_unsigned(max(-x, 0))

def signed_mul(x, y, normalize=True):
    xp, xn = x; yp, yn = y
    pp, nn = tensor(xp, yp), tensor(xn, yn)
    pn, npair = tensor(xp, yn), tensor(xn, yp)
    pos = direct_sum(pp, nn)
    neg = direct_sum(pn, npair)
    if normalize: pos, neg = norm_pair(pos), norm_pair(neg)
    return pos, neg

def signed_value(x):
    return val(x[0]) - val(x[1])

def algebra_check(seed=7, n=1000):
    rng = np.random.default_rng(seed)
    max_mul = max_add = max_gauge = 0.0
    for _ in range(n):
        a, b = 10**rng.uniform(-8, 8, 2)
        A, B = enc_unsigned(a), enc_unsigned(b)
        max_mul = max(max_mul, abs(val(norm_pair(tensor(A,B))) - a*b) / max(1, abs(a*b)))
        max_add = max(max_add, abs(val(norm_pair(direct_sum(A,B))) - (a+b)) / max(1, abs(a+b)))
        c = 10**rng.uniform(-10, 10)
        max_gauge = max(max_gauge, abs(val(A*c)-a) / max(1, abs(a)))
    max_signed = 0.0
    for _ in range(n):
        x, y = rng.uniform(-10, 10, 2)
        got = signed_value(signed_mul(enc_signed(x), enc_signed(y)))
        max_signed = max(max_signed, abs(got-x*y)/max(1,abs(x*y)))
    return dict(max_unsigned_mul_relerr=max_mul, max_unsigned_add_relerr=max_add,
                max_signed_mul_relerr=max_signed, max_gauge_relerr=max_gauge)

def product_stress(dtype=np.float32, depth=200, factor=1.25):
    # Positive products: normalized pair rails are bounded, but multiplication itself
    # still performs the same potentially overflowing products as the scalar control.
    f = np.array(factor, dtype=dtype)
    scalar = np.array(1.0, dtype=dtype)
    pair = np.array([1.0, factor], dtype=dtype)
    finite_scalar = True; finite_pair = True; max_rail = 0.0
    for _ in range(depth):
        scalar = scalar * f
        pair = pair * np.array([1.0, factor], dtype=dtype)
        s = np.max(np.abs(pair))
        if np.isfinite(s) and s > 0: pair = pair / s
        finite_scalar &= bool(np.isfinite(scalar))
        finite_pair &= bool(np.all(np.isfinite(pair)))
        max_rail = max(max_rail, float(np.max(np.abs(pair))) if np.all(np.isfinite(pair)) else math.inf)
    return dict(dtype=str(dtype), depth=depth, factor=factor,
                scalar_finite=finite_scalar, pair_finite=finite_pair,
                scalar_value=float(scalar) if np.isfinite(scalar) else None,
                pair_rails=pair.tolist() if np.all(np.isfinite(pair)) else None,
                decoded_pair=float(pair[1]/pair[0]) if pair[0] != 0 and np.all(np.isfinite(pair)) else None,
                max_normalized_rail=max_rail)

def cancellation_check():
    # A signed representation with both sign rails active has cancellation; compare
    # decoded error after repeated multiplication against ordinary scalar arithmetic.
    vals = [1.0 + 1e-5, 1.0, 1.0 - 1e-5]
    x = enc_signed(vals[0]); y = enc_signed(vals[1]); z = enc_signed(vals[2])
    pair = signed_mul(signed_mul(x,y),z)
    return dict(expected=math.prod(vals), decoded=signed_value(pair),
                abs_error=abs(signed_value(pair)-math.prod(vals)))

def main():
    np.seterr(over='ignore', under='ignore', invalid='ignore', divide='ignore')
    out = dict(algebra=algebra_check(), cancellation=cancellation_check(),
               stress=[product_stress(np.float32, 200, 1.25),
                       product_stress(np.float32, 100, 1.1),
                       product_stress(np.float64, 1000, 1.1)])
    print(json.dumps(out, indent=2, allow_nan=False))

if __name__ == '__main__': main()