g_math 0.6.8

#!/usr/bin/env python3
"""References for the compute-tier rework of the projective module, the
vector cross product and normalization, the tensor decompositions and the
pseudoinverse (tests/data/projective_tensor_compute_tier_refs.rs).

Inputs are dyadic with at most 8 fraction bits and small magnitude, so they
are stored exactly at every gated split (8 to 24 fraction bits) and on every
wider profile. References are exact rationals (fractions.Fraction), written
exactly when they terminate in decimal and otherwise through mpmath at 120
significant digits; the test parses each with the exact literal parser, which
rounds it to nearest at the build's split, so one table serves every build.

Usage: python3 scripts/generate_projective_tensor_compute_tier_refs.py
"""
import random
from fractions import Fraction as Fr

from mpmath import mp, mpf, sqrt

mp.dps = 140
SEED = 20260925
DIGITS = 120


def dy(rng, lo, hi, bits=8):
    """A dyadic k / 2^bits in [lo, hi)."""
    return Fr(rng.randrange(int(lo * (1 << bits)), int(hi * (1 << bits))), 1 << bits)


def dyn(rng, lo, hi, bits=8):
    """A nonzero dyadic in [lo, hi)."""
    while True:
        v = dy(rng, lo, hi, bits)
        if v != 0:
            return v


def lit(v):
    """A literal for a Fraction (exact when it terminates) or an mpf."""
    if isinstance(v, Fr):
        sign = "-" if v < 0 else ""
        a = abs(v)
        d = 0
        while (a * 10 ** d).denominator != 1:
            d += 1
            if d > 60:
                return sign + lit(mpf(a.numerator) / a.denominator)
        n = int(a * 10 ** d)
        if d == 0:
            return f"{sign}{n}"
        s = str(n).rjust(d + 1, "0")
        return f"{sign}{s[:-d]}.{s[-d:]}"
    return mp.nstr(v, DIGITS, strip_zeros=False, min_fixed=-10**9, max_fixed=10**9)


def sl(v):
    return "&[" + ", ".join(f'"{lit(x)}"' for x in v) + "]"


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


def cr1(a, b, c, d):
    return (a - c) * (b - d) / ((a - d) * (b - c))


def main():
    rng = random.Random(SEED)
    out = [
        "// GENERATED by scripts/generate_projective_tensor_compute_tier_refs.py - do not edit.",
        "// Dyadic inputs (<= 8 fraction bits); references exact or mpmath at 120",
        "// digits, parsed at the build's split (nearest, ties toward +inf).",
        "",
    ]

    # ---- cross_ratio_1d -------------------------------------------------
    rows = []
    while len(rows) < 12:
        a, b, c, d = (dy(rng, -4, 4) for _ in range(4))
        if len({a, b, c, d}) < 4:
            continue
        r = cr1(a, b, c, d)
        if abs(r) > 32:
            continue
        rows.append((a, b, c, d, r))
    # near-coincident points: small differences, bounded ratio
    while len(rows) < 24:
        a = dy(rng, -3, 3)
        c = a + Fr(rng.choice([-3, -2, -1, 1, 2, 3]), 256)
        d = a + Fr(rng.choice([-5, -4, -3, 3, 4, 5]), 256)
        b = c + Fr(rng.choice([-7, -3, -1, 1, 2, 5]), 256)
        if len({a, b, c, d}) < 4:
            continue
        r = cr1(a, b, c, d)
        if abs(r) > 32:
            continue
        rows.append((a, b, c, d, r))
    # (a - d)(b - c) = 2^-16: rounds to zero in storage at 8 fraction bits
    a = Fr(0)
    d = Fr(1, 256)
    c = Fr(2, 256)
    b = Fr(3, 256)
    rows.append((a, b, c, d, cr1(a, b, c, d)))
    out.append("/// (a, b, c, d, (a-c)(b-d) / ((a-d)(b-c)))")
    out.append("pub const CROSS_RATIO_1D: &[(&str, &str, &str, &str, &str)] = &[")
    for a, b, c, d, r in rows:
        out.append(f'    ("{lit(a)}", "{lit(b)}", "{lit(c)}", "{lit(d)}", "{lit(r)}"),')
    out.append("];\n")

    # ---- cross_ratio (collinear points in R^2 / R^3) ----------------------
    rows = []
    while len(rows) < 12:
        n = rng.choice([2, 3])
        a = [dy(rng, -2, 2, 4) for _ in range(n)]
        dirv = [dy(rng, -2, 2, 2) for _ in range(n)]
        if all(x == 0 for x in dirv):
            continue
        tc = Fr(rng.randrange(-12, 12), 4)
        td = Fr(rng.randrange(-12, 12), 4)
        if len({Fr(0), Fr(1), tc, td}) < 4:
            continue
        r = cr1(Fr(0), Fr(1), tc, td)
        if abs(r) > 32:
            continue
        b = [a[i] + dirv[i] for i in range(n)]
        c = [a[i] + tc * dirv[i] for i in range(n)]
        d = [a[i] + td * dirv[i] for i in range(n)]
        rows.append((a, b, c, d, r))
    # a short direction: |b - a|^2 = 2^-15 rounds to zero at 8 fraction bits
    a = [Fr(1, 2), Fr(-1, 4)]
    dirv = [Fr(1, 256), Fr(-1, 256)]
    for tc, td in [(Fr(2), Fr(-3)), (Fr(-1), Fr(3))]:
        b = [a[i] + dirv[i] for i in range(2)]
        c = [a[i] + tc * dirv[i] for i in range(2)]
        d = [a[i] + td * dirv[i] for i in range(2)]
        rows.append((a, b, c, d, cr1(Fr(0), Fr(1), tc, td)))
    out.append("/// (a, b, c, d, cross ratio) for collinear points")
    out.append("pub const CROSS_RATIO_ND: &[(&[&str], &[&str], &[&str], &[&str], &str)] = &[")
    for a, b, c, d, r in rows:
        out.append(f'    ({sl(a)}, {sl(b)}, {sl(c)}, {sl(d)}, "{lit(r)}"),')
    out.append("];\n")

    # ---- projective_transform ---------------------------------------------
    rows = []
    while len(rows) < 16:
        h = [dy(rng, -2, 2) for _ in range(9)]
        p = [dy(rng, -2, 2) for _ in range(2)]
        hp = p + [Fr(1)]
        y = [sum(h[3 * i + j] * hp[j] for j in range(3)) for i in range(3)]
        w = y[2]
        small = len(rows) >= 12
        if w == 0 or (abs(w) < Fr(1, 8)) != small or abs(w) < Fr(1, 64):
            continue
        x = [y[0] / w, y[1] / w]
        if max(abs(t) for t in x) > 64:
            continue
        rows.append((h, p, x))
    out.append("/// (H row-major 3x3, point, H applied and dehomogenized)")
    out.append("pub const PROJECTIVE: &[(&[&str], &[&str], &[&str])] = &[")
    for h, p, x in rows:
        out.append(f"    ({sl(h)}, {sl(p)}, {sl(x)}),")
    out.append("];\n")

    # ---- stereo_unproject -----------------------------------------------
    rows = []
    while len(rows) < 10:
        n = rng.choice([2, 3])
        x = [dy(rng, -4, 4) for _ in range(n)]
        s = sum(t * t for t in x)
        p = [2 * t / (1 + s) for t in x] + [(s - 1) / (s + 1)]
        rows.append((x, p))
    out.append("/// (x, inverse stereographic image)")
    out.append("pub const STEREO_UNPROJECT: &[(&[&str], &[&str])] = &[")
    for x, p in rows:
        out.append(f"    ({sl(x)}, {sl(p)}),")
    out.append("];\n")

    # ---- Moebius (real) -------------------------------------------------
    rows = []
    while len(rows) < 12:
        a, b, c, d, x = (dy(rng, -4, 4) for _ in range(5))
        den = c * x + d
        if abs(den) < Fr(1, 16):
            continue
        r = (a * x + b) / den
        if abs(r) > 32:
            continue
        rows.append(([a, b, c, d, x], r))
    out.append("/// ([a, b, c, d, x], (ax + b) / (cx + d))")
    out.append("pub const MOEBIUS_APPLY: &[(&[&str], &str)] = &[")
    for v, r in rows:
        out.append(f'    ({sl(v)}, "{lit(r)}"),')
    out.append("];\n")

    rows = []
    for _ in range(10):
        m1 = [dy(rng, -4, 4) for _ in range(4)]
        m2 = [dy(rng, -4, 4) for _ in range(4)]
        a, b, c, d = m1
        e, f, g, h = m2
        comp = [a * e + b * g, a * f + b * h, c * e + d * g, c * f + d * h]
        rows.append((m1, m2, comp, a * d - b * c))
    out.append("/// (m1, m2, m1 * m2, det m1)")
    out.append("pub const MOEBIUS_COMPOSE: &[(&[&str], &[&str], &[&str], &str)] = &[")
    for m1, m2, comp, det in rows:
        out.append(f'    ({sl(m1)}, {sl(m2)}, {sl(comp)}, "{lit(det)}"),')
    out.append("];\n")

    # ---- Moebius (complex) ----------------------------------------------
    def cmul(p, q):
        return (p[0] * q[0] - p[1] * q[1], p[0] * q[1] + p[1] * q[0])

    def cadd(p, q):
        return (p[0] + q[0], p[1] + q[1])

    rows = []
    while len(rows) < 12:
        v = [dy(rng, -2, 2) for _ in range(10)]
        a, b, c, d, z = [(v[2 * i], v[2 * i + 1]) for i in range(5)]
        num = cadd(cmul(a, z), b)
        den = cadd(cmul(c, z), d)
        nn = den[0] ** 2 + den[1] ** 2
        if nn < Fr(1, 16):
            continue
        r = ((num[0] * den[0] + num[1] * den[1]) / nn, (num[1] * den[0] - num[0] * den[1]) / nn)
        if max(abs(r[0]), abs(r[1])) > 32:
            continue
        rows.append((v, r))
    out.append("/// ([a, b, c, d, z] as (re, im) pairs, (az + b) / (cz + d))")
    out.append("pub const MOEBIUS_COMPLEX_APPLY: &[(&[&str], &[&str])] = &[")
    for v, r in rows:
        out.append(f"    ({sl(v)}, {sl(r)}),")
    out.append("];\n")

    rows = []
    for _ in range(8):
        v1 = [dy(rng, -2, 2) for _ in range(8)]
        v2 = [dy(rng, -2, 2) for _ in range(8)]
        a, b, c, d = [(v1[2 * i], v1[2 * i + 1]) for i in range(4)]
        e, f, g, h = [(v2[2 * i], v2[2 * i + 1]) for i in range(4)]
        comp = [cadd(cmul(a, e), cmul(b, g)), cadd(cmul(a, f), cmul(b, h)),
                cadd(cmul(c, e), cmul(d, g)), cadd(cmul(c, f), cmul(d, h))]
        rows.append((v1, v2, [t for p in comp for t in p]))
    out.append("/// (m1, m2, m1 o m2) as (re, im) pairs of a, b, c, d")
    out.append("pub const MOEBIUS_COMPLEX_COMPOSE: &[(&[&str], &[&str], &[&str])] = &[")
    for v1, v2, comp in rows:
        out.append(f"    ({sl(v1)}, {sl(v2)}, {sl(comp)}),")
    out.append("];\n")

    # ---- cross product and normalization ---------------------------------
    rows = []
    for _ in range(12):
        u = [dy(rng, -4, 4) for _ in range(3)]
        v = [dy(rng, -4, 4) for _ in range(3)]
        w = [u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2], u[0] * v[1] - u[1] * v[0]]
        rows.append((u, v, w))
    out.append("/// (u, v, u x v)")
    out.append("pub const CROSS: &[(&[&str], &[&str], &[&str])] = &[")
    for u, v, w in rows:
        out.append(f"    ({sl(u)}, {sl(v)}, {sl(w)}),")
    out.append("];\n")

    rows = []
    while len(rows) < 12:
        n = rng.randrange(2, 6)
        v = [dy(rng, -4, 4) for _ in range(n)]
        if all(t == 0 for t in v):
            continue
        rows.append(v)
    rows.append([Fr(1, 256), Fr(-3, 256), Fr(2, 256)])
    rows.append([Fr(5, 256), Fr(1, 256)])
    out.append("/// (v, v / |v|)")
    out.append("pub const NORMALIZE: &[(&[&str], &[&str])] = &[")
    for v in rows:
        nv = sqrt(sum(m_of(t) ** 2 for t in v))
        out.append(f"    ({sl(v)}, {sl([m_of(t) / nv for t in v])}),")
    out.append("];\n")

    # ---- TruncatedSVD::reconstruct ----------------------------------------
    rows = []
    for _ in range(8):
        m, k, n = rng.choice([(3, 2, 3), (4, 2, 3), (2, 2, 4), (3, 3, 3)])
        u = [dy(rng, -2, 2) for _ in range(m * k)]
        s = [dy(rng, 0, 4) for _ in range(k)]
        vt = [dy(rng, -2, 2) for _ in range(k * n)]
        r = [sum(u[i * k + t] * s[t] * vt[t * n + j] for t in range(k)) for i in range(m) for j in range(n)]
        rows.append((m, k, n, u, s, vt, r))
    out.append("/// (m, k, n, U m x k, sigma, Vt k x n, U diag(sigma) Vt)")
    out.append("pub const TSVD_RECONSTRUCT: &[(usize, usize, usize, &[&str], &[&str], &[&str], &[&str])] = &[")
    for m, k, n, u, s, vt, r in rows:
        out.append(f"    ({m}, {k}, {n}, {sl(u)}, {sl(s)}, {sl(vt)}, {sl(r)}),")
    out.append("];\n")

    # ---- Tucker reconstruct: core 2x2x2, factors 3x2 ----------------------
    rows = []
    for _ in range(6):
        g = [dy(rng, -2, 2) for _ in range(8)]
        fs = [[dy(rng, -1, 1) for _ in range(6)] for _ in range(3)]
        t = []
        for i in range(3):
            for j in range(3):
                for l in range(3):
                    t.append(sum(g[4 * p + 2 * q + r] * fs[0][2 * i + p] * fs[1][2 * j + q] * fs[2][2 * l + r]
                                 for p in range(2) for q in range(2) for r in range(2)))
        rows.append((g, fs, t))
    out.append("/// (core 2x2x2, U0 3x2, U1 3x2, U2 3x2, reconstructed 3x3x3)")
    out.append("pub const TUCKER_RECONSTRUCT: &[(&[&str], &[&str], &[&str], &[&str], &[&str])] = &[")
    for g, fs, t in rows:
        out.append(f"    ({sl(g)}, {sl(fs[0])}, {sl(fs[1])}, {sl(fs[2])}, {sl(t)}),")
    out.append("];\n")

    # ---- CP reconstruct: rank 2, factors 3x2 --------------------------------
    rows = []
    for _ in range(6):
        w = [dy(rng, 0, 4) for _ in range(2)]
        fs = [[dy(rng, -1, 1) for _ in range(6)] for _ in range(3)]
        t = []
        for i in range(3):
            for j in range(3):
                for l in range(3):
                    t.append(sum(w[r] * fs[0][2 * i + r] * fs[1][2 * j + r] * fs[2][2 * l + r] for r in range(2)))
        rows.append((w, fs, t))
    out.append("/// (weights, A0 3x2, A1 3x2, A2 3x2, reconstructed 3x3x3)")
    out.append("pub const CP_RECONSTRUCT: &[(&[&str], &[&str], &[&str], &[&str], &[&str])] = &[")
    for w, fs, t in rows:
        out.append(f"    ({sl(w)}, {sl(fs[0])}, {sl(fs[1])}, {sl(fs[2])}, {sl(t)}),")
    out.append("];\n")

    # ---- CP-ALS on exactly rank-1 tensors ---------------------------------
    # T = w a o b o c. In exact arithmetic one ALS sweep from any start not
    # orthogonal to the factors returns directions a, b, c and the weight
    # |w| |a| |b| |c|, independent of the start: the SVD initialization
    # cancels, so the reference is deterministic.
    rows = []
    while len(rows) < 6:
        w = dyn(rng, Fr(1, 2), 2, 2)
        vs = [[dy(rng, -1, 1, 2) for _ in range(3)] for _ in range(3)]
        if any(all(x == 0 for x in v) for v in vs):
            continue
        t = [w * vs[0][i] * vs[1][j] * vs[2][l] for i in range(3) for j in range(3) for l in range(3)]
        norms = [sqrt(sum(m_of(x) ** 2 for x in v)) for v in vs]
        weight = abs(m_of(w)) * norms[0] * norms[1] * norms[2]
        dirs = [[abs(m_of(x)) / norms[k] for x in vs[k]] for k in range(3)]
        rows.append((t, weight, dirs))
    out.append("/// (rank-1 tensor 3x3x3, weight, |a|/||a||, |b|/||b||, |c|/||c||)")
    out.append("pub const CP_RANK1: &[(&[&str], &str, &[&str], &[&str], &[&str])] = &[")
    for t, weight, dirs in rows:
        out.append(f'    ({sl(t)}, "{lit(weight)}", {sl(dirs[0])}, {sl(dirs[1])}, {sl(dirs[2])}),')
    out.append("];\n")

    # ---- pseudoinverse ------------------------------------------------------
    # A = U diag(s) V^T with U, V signed 4x4 Hadamard / 2 (orthogonal, dyadic)
    # and dyadic singular values, so A is dyadic and A+ = V diag(1/s) U^T.
    had = [[1, 1, 1, 1], [1, -1, 1, -1], [1, 1, -1, -1], [1, -1, -1, 1]]
    rows = []
    for case in range(6):
        su = [rng.choice([1, -1]) for _ in range(4)]
        sv = [rng.choice([1, -1]) for _ in range(4)]
        pu = list(range(4))
        pv = list(range(4))
        rng.shuffle(pu)
        rng.shuffle(pv)
        uu = [[Fr(had[i][pu[j]] * su[j], 2) for j in range(4)] for i in range(4)]
        vv = [[Fr(had[i][pv[j]] * sv[j], 2) for j in range(4)] for i in range(4)]
        if case < 2:
            s = [Fr(3), Fr(0), Fr(0), Fr(0)]  # rank 1
        else:
            s = sorted([Fr(rng.randrange(8, 64), 16) for _ in range(4)], reverse=True)
            if len(set(s)) < 4:
                s = [Fr(3), Fr(2), Fr(3, 2), Fr(3, 4)]
        a = [sum(uu[i][t] * s[t] * vv[j][t] for t in range(4)) for i in range(4) for j in range(4)]
        p = [sum(vv[i][t] * (1 / s[t]) * uu[j][t] for t in range(4) if s[t] != 0) for i in range(4) for j in range(4)]
        rows.append((4, 4, a, p))
    # signed permutations times a diagonal: the SVD factors are exact, so the
    # pseudoinverse's own arithmetic is all that is measured
    for _ in range(4):
        perm = list(range(4))
        rng.shuffle(perm)
        s = [Fr(rng.choice([3, 5, 7, 9, 11, 13]), rng.choice([1, 2, 4])) * rng.choice([1, -1]) for _ in range(4)]
        a = [[Fr(0)] * 4 for _ in range(4)]
        for i in range(4):
            a[i][perm[i]] = s[i]
        # A = P D with P a permutation: A+ = A^-1 (entries 1/s at the transposed position)
        p = [[Fr(0)] * 4 for _ in range(4)]
        for i in range(4):
            p[perm[i]][i] = 1 / s[i]
        rows.append((4, 4, [x for r in a for x in r], [x for r in p for x in r]))
    out.append("/// (m, n, A row-major, pseudoinverse n x m row-major); the last four")
    out.append("/// are signed permutations of a diagonal (exact SVD factors)")
    out.append("pub const PINV: &[(usize, usize, &[&str], &[&str])] = &[")
    for m, n, a, p in rows:
        out.append(f"    ({m}, {n}, {sl(a)}, {sl(p)}),")
    out.append("];\n")

    # ---- condition_number_1 -------------------------------------------------
    def inv(mat):
        n = len(mat)
        a = [row[:] + [Fr(int(i == j)) for j in range(n)] for i, row in enumerate(mat)]
        for k in range(n):
            piv = next(i for i in range(k, n) if a[i][k] != 0)
            a[k], a[piv] = a[piv], a[k]
            pk = a[k][k]
            a[k] = [x / pk for x in a[k]]
            for i in range(n):
                if i != k and a[i][k] != 0:
                    f = a[i][k]
                    a[i] = [x - f * y for x, y in zip(a[i], a[k])]
        return [row[n:] for row in a]

    def norm1(mat):
        n = len(mat)
        return max(sum(abs(mat[i][j]) for i in range(n)) for j in range(n))

    rows = []
    while len(rows) < 6:
        n = 3
        mat = [[dy(rng, -2, 2, 4) for _ in range(n)] for _ in range(n)]
        for i in range(n):
            mat[i][i] += 3
        try:
            mi = inv(mat)
        except StopIteration:
            continue
        rows.append((n, [x for r in mat for x in r], norm1(mat) * norm1(mi)))
    out.append("/// (n, A row-major, ||A||_1 ||A^-1||_1)")
    out.append("pub const COND1: &[(usize, &[&str], &str)] = &[")
    for n, a, c in rows:
        out.append(f'    ({n}, {sl(a)}, "{lit(c)}"),')
    out.append("];")

    with open("tests/data/projective_tensor_compute_tier_refs.rs", "w") as fh:
        fh.write("\n".join(out) + "\n")


if __name__ == "__main__":
    main()