#!/usr/bin/env python3
"""Precompute mpmath/MPFR O3 tables for math calculators that share the arccos issue.

Re-run: python3 scripts/lib/cvp/oracles/generate-math-o3.py
Public copies: /developers/cvp/reproduce/generate-math-o3.py and math-o3-tables.json
"""
from __future__ import annotations

import json
import math
import sys
from pathlib import Path

import mpmath as mp

DPS = 80
SEED = 20260907
GENERATOR_ID = "math-mpmath-o3"
GENERATOR_VERSION = "1.0.0"


def lcg(seed: int):
    x = seed & 0xFFFFFFFF

    def nxt() -> float:
        nonlocal x
        x = (1664525 * x + 1013904223) & 0xFFFFFFFF
        return x / 2**32

    return nxt


def vec(vid: str, inputs: dict, y: mp.mpf, *, kind: str, extra: dict | None = None) -> dict:
    row = {
        "id": vid,
        "kind": kind,
        "inputs": inputs,
        "y_f64": float(y),
        "y_decimal": mp.nstr(y, 40, strip_zeros=False),
    }
    if extra:
        row.update(extra)
    return row


def family_meta(name: str, notes: str, vectors: list[dict]) -> dict:
    return {
        "family": name,
        "generator_id": GENERATOR_ID,
        "generator_version": GENERATOR_VERSION,
        "seed": str(SEED),
        "mpmath_dps": DPS,
        "precision_bits": int(DPS * math.log2(10)),
        "library": f"mpmath {mp.__version__}",
        "notes": notes,
        "vectors": vectors,
    }


def invtrig(fn: str) -> list[dict]:
    mp_fn = {"asin": mp.asin, "atan": mp.atan}[fn]
    named = []
    if fn == "asin":
        xs = [
            ("neg-one", -1.0, True),
            ("one", 1.0, True),
            ("zero", 0.0, True),
            ("half", 0.5, True),
            ("neg-half", -0.5, True),
            ("sqrt2-2", math.sqrt(0.5), True),
            ("neg-sqrt2-2", -math.sqrt(0.5), True),
            ("no-snap", 0.5000000009, False),
            ("near-one", 1.0 - 2.0**-40, False),
            ("near-neg-one", -1.0 + 2.0**-40, False),
            ("quarter", 0.25, False),
        ]
    else:
        xs = [
            ("zero", 0.0, True),
            ("one", 1.0, True),
            ("neg-one", -1.0, True),
            ("sqrt3", math.sqrt(3), True),
            ("inv-sqrt3", 1 / math.sqrt(3), True),
            ("no-snap", 1.0000000009, False),
            ("ten", 10.0, False),
            ("neg-ten", -10.0, False),
            ("quarter", 0.25, False),
            ("tiny", 1e-20, False),
            ("1e10", 1e10, False),
            ("1e20", 1e20, False),
            ("max", sys.float_info.max, False),
        ]
    for name, x, exact in xs:
        rad = mp_fn(mp.mpf(x))
        deg = rad * mp.mpf(180) / mp.pi
        named.append(
            vec(
                f"o3-{fn}-{name}",
                {"x": x},
                deg,
                kind="algebraic" if exact else "interior",
                extra={
                    "exact_algebraic": exact,
                    "rad_f64": float(rad),
                    "rad_decimal": mp.nstr(rad, 40, strip_zeros=False),
                },
            )
        )
    rnd = lcg(SEED + (1 if fn == "asin" else 2))
    seen = {v["inputs"]["x"] for v in named}
    i = 0
    while i < 8:
        u = rnd()
        x = (2.0 * u - 1.0) if fn == "asin" else math.tan((u - 0.5) * 1.2)
        if fn == "asin" and abs(x) > 0.999:
            continue
        if x in seen or not math.isfinite(x):
            continue
        rad = mp_fn(mp.mpf(x))
        deg = rad * mp.mpf(180) / mp.pi
        named.append(
            vec(
                f"o3-{fn}-rand-{i:02d}",
                {"x": x},
                deg,
                kind="random-interior",
                extra={
                    "exact_algebraic": False,
                    "rad_f64": float(rad),
                    "rad_decimal": mp.nstr(rad, 40, strip_zeros=False),
                },
            )
        )
        seen.add(x)
        i += 1
    return named


def exact_trig(fn: str, deg: float) -> mp.mpf | None:
    """Algebraic values the IUT snaps — do not go through π, or cos(90°) becomes ~1e-81."""
    half = mp.mpf("0.5")
    s2 = mp.sqrt(2) / 2
    s3 = mp.sqrt(3) / 2
    t3 = mp.sqrt(3)
    table = {
        "sin": {0: mp.mpf(0), 30: half, 45: s2, 60: s3, 90: mp.mpf(1), 180: mp.mpf(0), -30: -half},
        "cos": {0: mp.mpf(1), 30: s3, 45: s2, 60: half, 90: mp.mpf(0), 180: mp.mpf(-1), -30: s3},
        "tan": {0: mp.mpf(0), 30: 1 / t3, 45: mp.mpf(1), 60: t3, -30: -1 / t3},
    }
    return table.get(fn, {}).get(deg)


def trig(fn: str) -> list[dict]:
    mp_fn = {"sin": mp.sin, "cos": mp.cos, "tan": mp.tan}[fn]
    angles = [0, 30, 45, 60]
    if fn != "tan":
        angles += [90, 180]
    angles += [33.3, 12.5, -30, 15]
    out = []
    for deg in angles:
        if fn == "tan" and abs((deg - 90) % 180) < 1e-12:
            continue
        y = exact_trig(fn, deg)
        if y is None:
            rad = mp.mpf(deg) * mp.pi / mp.mpf(180)
            y = mp_fn(rad)
        out.append(vec(f"o3-{fn}-{deg}", {"theta": float(deg), "unit": "deg"}, y, kind="angle"))
    rnd = lcg(SEED + {"sin": 11, "cos": 12, "tan": 13}[fn])
    for i in range(8):
        deg = (rnd() * 100) - 50  # stay inside (-50, 50) — tan ULP stays small vs Node
        if fn == "tan" and abs((deg - 90) % 180) < 1:
            continue
        rad = mp.mpf(deg) * mp.pi / mp.mpf(180)
        y = mp_fn(rad)
        out.append(vec(f"o3-{fn}-rand-{i:02d}", {"theta": float(deg), "unit": "deg"}, y, kind="random-interior"))
    return out


def ln_family() -> list[dict]:
    xs = [1.0, math.e, 10.0, 0.5, 2.0, 100.0, 1e-6, 1e6]
    out = []
    for x in xs:
        y = mp.mpf(0) if x == 1.0 else mp.mpf(1) if x == math.e else mp.ln(mp.mpf(x))
        tag = str(x).replace(".", "p").replace("-", "m")
        out.append(vec(f"o3-ln-{tag}", {"x": x}, y, kind="interior"))
    rnd = lcg(SEED + 20)
    for i in range(8):
        x = 10 ** (rnd() * 6 - 3)
        y = mp.ln(mp.mpf(x))
        out.append(vec(f"o3-ln-rand-{i:02d}", {"x": float(x)}, y, kind="random-interior"))
    return out


def log10_family() -> list[dict]:
    xs = [1.0, 10.0, 100.0, 0.1, 2.0, 1000.0]
    exact = {1.0: 0, 10.0: 1, 100.0: 2, 0.1: -1, 1000.0: 3}
    out = []
    for x in xs:
        y = mp.mpf(exact[x]) if x in exact else mp.log10(mp.mpf(x))
        tag = str(x).replace(".", "p")
        out.append(vec(f"o3-log10-{tag}", {"x": x, "base": "10"}, y, kind="interior"))
    rnd = lcg(SEED + 21)
    for i in range(8):
        x = 10 ** (rnd() * 6 - 2)
        y = mp.log10(mp.mpf(x))
        out.append(vec(f"o3-log10-rand-{i:02d}", {"x": float(x), "base": "10"}, y, kind="random-interior"))
    return out


def sqrt_family() -> list[dict]:
    half = mp.mpf("0.5")
    xs = [0.0, 1.0, 4.0, 9.0, 144.0, 2.0, 0.25, 1e-8, 1e8]
    exact = {0.0: 0, 1.0: 1, 4.0: 2, 9.0: 3, 144.0: 12, 0.25: half}
    out = []
    for x in xs:
        y = mp.mpf(exact[x]) if x in exact else mp.sqrt(mp.mpf(x))
        tag = str(x).replace(".", "p")
        out.append(vec(f"o3-sqrt-{tag}", {"x": x}, y, kind="interior"))
    rnd = lcg(SEED + 22)
    for i in range(8):
        x = rnd() * 1e4
        y = mp.sqrt(mp.mpf(x))
        out.append(vec(f"o3-sqrt-rand-{i:02d}", {"x": float(x)}, y, kind="random-interior"))
    return out


def exponent_family() -> list[dict]:
    pairs = [(2, 10), (10, 2), (5, 0), (9, 0.5), (2, -3), (math.e, 1)]
    out = []
    for a, n in pairs:
        y = mp.e ** mp.mpf(n) if abs(a - math.e) < 1e-12 else mp.mpf(a) ** mp.mpf(n)
        out.append(vec(f"o3-exp-{a}-{n}", {"a": a, "n": n}, y, kind="interior"))
    return out


def antilog_family() -> list[dict]:
    """x = b^y with the same bases the IUT accepts as tokens.

    Integer powers prove basic consistency. Extra vectors challenge binary64:
    10^0.5, 2^(≈1/3), e^3, near-base-one, max-finite, and subnormal/tiny.
    """
    out = []
    for y in (0.0, 1.0, 2.0, 3.0, -1.0):
        val = mp.mpf(10) ** mp.mpf(y)
        out.append(vec(f"o3-antilog-10-{y}", {"y": float(y), "base": "10"}, val, kind="interior"))
    out.append(vec("o3-antilog-2-3", {"y": 3.0, "base": "2"}, mp.mpf(8), kind="interior"))
    out.append(vec("o3-antilog-e-1", {"y": 1.0, "base": "e"}, mp.e, kind="interior"))
    out.append(vec("o3-antilog-e-0", {"y": 0.0, "base": "e"}, mp.mpf(1), kind="interior"))
    out.append(
        vec(
            "o3-antilog-10-half",
            {"y": 0.5, "base": "10"},
            mp.mpf(10) ** mp.mpf("0.5"),
            kind="interior",
        )
    )
    y_cbrt = 0.3333333333333333  # IEEE-754 binary64 1/3, matching IUT input
    out.append(
        vec(
            "o3-antilog-2-cbrt",
            {"y": y_cbrt, "base": "2"},
            mp.mpf(2) ** mp.mpf(str(y_cbrt)),
            kind="interior",
        )
    )
    out.append(vec("o3-antilog-e-3", {"y": 3.0, "base": "e"}, mp.e ** 3, kind="interior"))
    near = 1.000000000000001
    out.append(
        vec(
            "o3-antilog-near-base-one",
            {"y": 2.0, "b": near},
            mp.mpf(str(near)) ** 2,
            kind="interior",
        )
    )
    out.append(
        vec(
            "o3-antilog-10-308",
            {"y": 308.0, "base": "10"},
            mp.mpf(10) ** 308,
            kind="boundary-finite",
        )
    )
    out.append(
        vec(
            "o3-antilog-10-neg323",
            {"y": -323.0, "base": "10"},
            mp.mpf(10) ** -323,
            kind="subnormal",
        )
    )
    out.append(
        vec(
            "o3-antilog-tiny-base",
            {"y": 2.0, "b": 1e-100},
            mp.mpf("1e-100") ** 2,
            kind="interior",
        )
    )
    return out


def root_family() -> list[dict]:
    pairs = [(8, 3), (16, 4), (100, 2), (27, 3), (32, 5), (1, 7)]
    out = []
    for a, n in pairs:
        y = mp.mpf(a) ** (mp.mpf(1) / n)
        out.append(vec(f"o3-root-{a}-{n}", {"a": a, "n": n}, y, kind="interior"))
    return out


def exp_growth_family() -> list[dict]:
    """Periodic x0*(1+r/100)**t, continuous x0*e**(k*t), and inverse t/x0/r."""
    out = []
    periodic = [
        (100.0, 5.0, 2.0),
        (50.0, 4.0, 90.0),
        (100.0, -10.0, 10.0),
        (100.0, 5.0, -10.0),
        (0.0, 4.0, 90.0),
        (16.0, 100.0, 1.0),
        (10000.0, 5.0, 11.0),
    ]
    for x0, r, t in periodic:
        y = mp.mpf(x0) * (1 + mp.mpf(r) / 100) ** mp.mpf(t)
        tag = f"{x0:g}-{r:g}-{t:g}".replace(".", "p").replace("-", "m")
        out.append(
            vec(
                f"o3-eg-{tag}",
                {"x0": float(x0), "ratePct": float(r), "t": float(t)},
                y,
                kind="forward-periodic",
            )
        )
    y_cont = mp.mpf(50) * mp.e ** (mp.mpf("0.04") * mp.mpf(10))
    out.append(
        vec(
            "o3-eg-cont-50-0p04-10",
            {"x0": 50.0, "ratePct": 0.04, "t": 10.0, "model": "continuous"},
            y_cont,
            kind="forward-continuous",
        )
    )
    y_cont_neg = mp.mpf(50) * mp.e ** (mp.mpf("-0.04") * mp.mpf(10))
    out.append(
        vec(
            "o3-eg-cont-50-m0p04-10",
            {"x0": 50.0, "ratePct": -0.04, "t": 10.0, "model": "continuous"},
            y_cont_neg,
            kind="forward-continuous",
        )
    )
    y_t = mp.log(mp.mpf(3)) / mp.log(mp.mpf("1.05"))
    out.append(
        vec(
            "o3-eg-solve-t",
            {"x0": 10000.0, "ratePct": 5.0, "xt": 30000.0},
            y_t,
            kind="inverse-time",
            extra={"engine": "exponential-growth-time"},
        )
    )
    out.append(
        vec(
            "o3-eg-solve-x0",
            {"ratePct": 100.0, "t": 1.0, "xt": 32.0},
            mp.mpf(16),
            kind="inverse-x0",
            extra={"engine": "exponential-growth-x0"},
        )
    )
    out.append(
        vec(
            "o3-eg-solve-rate",
            {"x0": 16.0, "t": 1.0, "xt": 32.0},
            mp.mpf(100),
            kind="inverse-rate",
            extra={"engine": "exponential-growth-rate"},
        )
    )
    y_half = mp.log(mp.mpf("0.5")) / mp.log(mp.mpf("0.5"))
    out.append(
        vec(
            "o3-eg-half-life-50pct",
            {"x0": 100.0, "ratePct": -50.0, "xt": 50.0},
            y_half,
            kind="inverse-time",
            extra={"engine": "exponential-growth-time"},
        )
    )
    return out


def json_int(n: int):
    max_safe = 2**53 - 1
    if -max_safe <= n <= max_safe:
        return n
    return str(n)


def divide_family() -> list[dict]:
    """Quotient binary64 vs mpmath, plus exact Python integer long division (independent of JS BigInt)."""
    out = []
    pairs = [
        ("1-3", 1.0, 3.0),
        ("100-7", 100.0, 7.0),
        ("1p005", 1.005, 1.0),
        ("tiny", 1.0, 1e-300),
        ("max-safe", float(2**53 - 1), 7.0),
    ]
    for tag, a, b in pairs:
        y = mp.mpf(a) / mp.mpf(b)
        inputs = {"a": a, "b": b, "mode": "quotient"}
        if tag == "1p005":
            inputs["decimals"] = 2
        out.append(vec(f"o3-div-{tag}", inputs, y, kind="quotient"))

    def long_vec(vid: str, a: int | str, b: int | str, *, kind: str) -> dict:
        aa = int(a)
        bb = int(b)
        q, r = divmod(aa, bb)
        max_safe = 2**53 - 1
        y = float(aa) / float(bb) if 0 <= aa <= max_safe and 0 < bb <= max_safe else None
        return {
            "id": vid,
            "kind": kind,
            "inputs": {"a": json_int(aa), "b": json_int(bb), "mode": "long"},
            "y_f64": y,
            "y_decimal": None if y is None else mp.nstr(mp.mpf(y), 40, strip_zeros=False),
            "quotient": json_int(q),
            "remainder": json_int(r),
        }

    out.append(long_vec("o3-div-long-100-7", 100, 7, kind="long-normal"))
    out.append(long_vec("o3-div-long-past-safe", "9007199254740993", "1", kind="long-large-integer"))
    out.append(long_vec("o3-div-long-huge", "9007199254740993123456789", "7", kind="long-large-integer"))
    out.append(long_vec("o3-div-long-48", "1" + "0" * 47, "3", kind="long-boundary"))
    return out


def main() -> None:
    mp.mp.dps = DPS
    tables = {
        "asin": family_meta("asin", "mpmath.asin → deg/rad, independent of Node Math.asin", invtrig("asin")),
        "atan": family_meta("atan", "mpmath.atan → deg/rad, independent of Node Math.atan", invtrig("atan")),
        "sin": family_meta("sin", "mpmath.sin of degree inputs, independent of Node Math.sin", trig("sin")),
        "cos": family_meta("cos", "mpmath.cos of degree inputs, independent of Node Math.cos", trig("cos")),
        "tan": family_meta("tan", "mpmath.tan of degree inputs away from poles", trig("tan")),
        "ln": family_meta("ln", "mpmath.ln, independent of Node Math.log", ln_family()),
        "log10": family_meta("log10", "mpmath.log10, independent of Node change-of-base", log10_family()),
        "sqrt": family_meta("sqrt", "mpmath.sqrt, independent of Node Math.sqrt", sqrt_family()),
        "exponent": family_meta("exponent", "mpmath a**n, independent of Node ** / Math.exp", exponent_family()),
        "antilog": family_meta("antilog", "mpmath b**y / e**y, independent of Node ** / Math.exp", antilog_family()),
        "root": family_meta("root", "mpmath a**(1/n), independent of Node Math.pow/cbrt", root_family()),
        "exp_growth": family_meta(
            "exp_growth",
            "mpmath periodic (1+r/100)**t / e**(kt) and inverse t/x0/r, independent of Node ** / Math.exp",
            exp_growth_family(),
        ),
        "divide": family_meta(
            "divide",
            "mpmath a/b rounded to nearest binary64, plus Python int divmod for long division — independent of Node / and JS BigInt",
            divide_family(),
        ),
    }
    out = Path(__file__).with_name("math-o3-tables.json")
    payload = {
        "generator_id": GENERATOR_ID,
        "generator_version": GENERATOR_VERSION,
        "seed": str(SEED),
        "families": tables,
    }
    out.write_text(json.dumps(payload, indent=2) + "\n", encoding="utf-8")
    counts = {k: len(v["vectors"]) for k, v in tables.items()}
    print(f"wrote {out} {counts}")


if __name__ == "__main__":
    main()