"""Generate test-only certified Gauss brackets from exact rational arithmetic.

Newton/Decimal produces candidates only. Exact Fraction signs prove n distinct
Legendre roots; rational interval recurrence encloses the associated weights.
No runtime solver dependency, reference modification, or libm ULP premise.
"""
from decimal import Decimal, localcontext
from fractions import Fraction as F
from pathlib import Path
import math
import re
import json


def add(a, b):
    return a[0] + b[0], a[1] + b[1]


def neg(a):
    return -a[1], -a[0]


def mul(a, b):
    v = [x * y for x in a for y in b]
    return min(v), max(v)


def div(a, b):
    assert b[0] > 0 or b[1] < 0
    return mul(a, (1 / b[1], 1 / b[0]))


def scale(a, b):
    return mul(a, (F(b), F(b)))


def legendre(n, x):
    p, q = (F(1), F(1)), x
    for k in range(2, n + 1):
        p, q = q, scale(add(scale(mul(x, q), 2 * k - 1), neg(scale(p, k - 1))), F(1, k))
    derivative = div(scale(add(mul(x, q), neg(p)), n), add(mul(x, x), (-F(1), -F(1))))
    return q, derivative


def bracket(value):
    return math.nextafter(float(value[0]), -math.inf), math.nextafter(float(value[1]), math.inf)


def rule(n):
    result = []
    for i in range(n):
        with localcontext() as context:
            context.prec = 90
            z = Decimal(str(math.cos(math.pi * (i + .75) / (n + .5))))
            for _ in range(30):
                p, q = Decimal(1), z
                for k in range(2, n + 1):
                    p, q = q, ((2 * k - 1) * z * q - (k - 1) * p) / k
                d = n * (z * q - p) / (z * z - 1)
                z -= q / d
        lo, hi = bracket((F(z), F(z)))
        pl = legendre(n, (F(lo), F(lo)))[0][0]
        ph = legendre(n, (F(hi), F(hi)))[0][0]
        assert pl * ph < 0
        x = F(lo), F(hi)
        derivative = legendre(n, x)[1]
        w = div((F(2), F(2)), mul(add((F(1), F(1)), neg(mul(x, x))), mul(derivative, derivative)))
        result.append((lo, hi, *bracket(w)))
    result.sort()
    assert all(a[1] < b[0] for a, b in zip(result, result[1:]))
    assert sum(F(r[2]) for r in result) < 2 < sum(F(r[3]) for r in result)
    return result


def slab_fixtures():
    # Approved Appendix A supplies the independently implemented Cholesky /
    # Jacobi spectral candidates; Appendix E supplies certified transcendentals.
    document = Path(__file__).resolve().parents[3] / "docs/shell-transient-heat-transfer/numerical-oracle-scratch.md"
    blocks = re.findall(r"```python\n(.*?)```", document.read_text(encoding="utf-8-sig"), re.S)
    spectral = {}
    exec(blocks[0].split('out={"lambda"')[0], spectral)
    interval_scope = {}
    exec(blocks[-1].split("\nfor t in [.001")[0], interval_scope)
    iv, ia, im, isc = (interval_scope[k] for k in ("interval", "add", "mul", "scale"))
    ex = interval_scope["exneg"]
    n = 4

    def exact_inverse(matrix):
        a = [list(row) + [F(i == j) for j in range(n)] for i, row in enumerate(matrix)]
        for k in range(n):
            pivot = next(i for i in range(k, n) if a[i][k])
            a[k], a[pivot] = a[pivot], a[k]
            divisor = a[k][k]
            a[k] = [x / divisor for x in a[k]]
            for i in range(n):
                if i != k:
                    factor = a[i][k]
                    a[i] = [x - factor * y for x, y in zip(a[i], a[k])]
        return [row[n:] for row in a]

    def exact_mm(a, b):
        return [[sum((x*y for x, y in zip(row, col)), F(0)) for col in zip(*b)] for row in a]

    def matrix_multiply(a, b):
        return [[sum_intervals([im(x, y) for x, y in zip(row, col)]) for col in zip(*b)] for row in a]

    def sum_intervals(values):
        total = iv(0)
        for v in values:
            total = ia(total, v)
        return total

    def apply(a, v):
        return [sum_intervals([im(x, y) for x, y in zip(row, v)]) for row in a]

    def power(a, count):
        result = [[iv(i == j) for j in range(n)] for i in range(n)]
        while count:
            if count & 1:
                result = matrix_multiply(result, a)
            count >>= 1
            if count:
                a = matrix_multiply(a, a)
        return result

    c = [[F(0) for _ in range(5)] for _ in range(5)]
    k = [[F(0) for _ in range(5)] for _ in range(5)]
    for offset in (0, 2):
        for i in range(3):
            for j in range(3):
                c[offset+i][offset+j] += F([[4, 2, -1], [2, 16, 2], [-1, 2, 4]][i][j], 60)
                k[offset+i][offset+j] += F(2, 3)*[[7, -8, 1], [-8, 16, -8], [1, -8, 7]][i][j]
    cf, kf = [r[1:] for r in c[1:]], [r[1:] for r in k[1:]]
    inverse = exact_inverse(cf)
    d = exact_mm(inverse, kf)
    y = [-1-sum(inverse[i][j]*c[j+1][0] for j in range(n)) for i in range(n)]

    def matrix_exponential(t):
        # An independent enclosure of exp(-Cff^-1 Kff t) certifies all spectral
        # eigenpair/coefficient errors together, rather than assuming libm ULPs.
        # For ||B||inf<=1/2, the 64-term exponential remainder is bounded by
        # 2*(1/2)^65/65!. Entrywise intervals preserve subsequent squarings.
        magnitude = max(sum(abs(x*t) for x in row) for row in d)
        shift = 0
        while magnitude > F(1, 2):
            magnitude /= 2
            shift += 1
        a = [[iv(-x*t/F(2)**shift) for x in row] for row in d]
        term = [[iv(i == j) for j in range(n)] for i in range(n)]
        result = [row[:] for row in term]
        for order in range(1, 65):
            term = [[isc(x, F(1, order)) for x in row] for row in matrix_multiply(term, a)]
            result = [[ia(x, z) for x, z in zip(row, terms)] for row, terms in zip(result, term)]
        remainder = 2*F(1, 2)**65/math.factorial(65)
        result = [[ia(x, iv(-remainder, remainder)) for x in row] for row in result]
        for _ in range(shift):
            result = matrix_multiply(result, result)
        return result

    def outputs(error):
        temperature = [iv(1)] + [ia(iv(1), z) for z in error]
        derivative = [iv(0)] + apply([[iv(-x) for x in row] for row in d], error)
        flux = []
        for offset, s in ((0, -1), (0, 0), (0, 1), (2, -1), (2, 0), (2, 1)):
            deriv = [F(4)*F(s-F(1, 2)), F(-8*s), F(4)*F(s+F(1, 2))]
            flux.append(sum_intervals([isc(temperature[offset+i], -deriv[i]) for i in range(3)]))
        energy = sum_intervals([isc(temperature[j], sum(c[i][j] for i in range(5))) for j in range(5)])
        reaction = sum_intervals([isc(temperature[j], k[0][j]) for j in range(5)] +
                                 [isc(derivative[j], c[0][j]) for j in range(5)])
        return temperature + flux + [energy, reaction]

    def spectral_candidate(t, count):
        amplitudes = []
        for eigenvalue, coefficient in zip(spectral["lam"], spectral["coef"]):
            decay = ex(iv(F(eigenvalue)*t)) if count is None else power_scalar(
                interval_scope["div"](iv(1), iv(1+F(eigenvalue)*t/count)), count)
            amplitudes.append(isc(decay, F(coefficient)))
        return outputs(apply([[iv(F(x)) for x in row] for row in spectral["W"]], amplitudes))

    def power_scalar(x, count):
        result = iv(1)
        while count:
            if count & 1:
                result = im(result, x)
            count >>= 1
            if count:
                x = im(x, x)
        return result

    table = [[.141, 10.7, .0304, 10.9], [.0153, 1.28, .00468, .824],
             [.00598, .253, .00130, .00632], [.00139, .136, .000666, .00953],
             [.000142, .0284, .0000945, .000526], [.0000414, .00816, .0000159, .0000479]]
    fixtures, metrics = [], []
    for time, space in zip((.001, .01, .05, .1, .5, 1.), table):
        t = F(time)
        temperatures, flux, energy, _, _ = interval_scope["series_interval"](time)
        continuum = temperatures + [flux[i] for i in (0, 1, 2, 2, 3, 4)] + [energy, flux[0]]
        solutions = []
        spectral_enclosure = 0.
        for count in (None, 4096, 8192):
            if count is None:
                operator = matrix_exponential(t)
            else:
                a = [[cf[i][j]+t/count*kf[i][j] for j in range(n)] for i in range(n)]
                operator = power([[iv(x) for x in row] for row in exact_mm(exact_inverse(a), cf)], count)
            exact = outputs(apply(operator, [iv(x) for x in y]))
            candidate = spectral_candidate(t, count)
            enclosed = [(min(a[0], b[0]), max(a[1], b[1])) for a, b in zip(exact, candidate)]
            spectral_enclosure = max(spectral_enclosure, max(float(hi-lo) for lo, hi in enclosed))
            solutions.append(enclosed)
        fixtures.append((time, space, continuum, *solutions))
        metrics.append({"time": time, "spectral_and_matrix_function_enclosure": spectral_enclosure})
    print(json.dumps({"eigenvalues": spectral["lam"], "slab_oracle_metrics": metrics}, indent=2))
    return fixtures


if __name__ == "__main__":
    output = ["// Generated by thermal_oracle_generator.py; exact root/weight certificates.",
              "#ifndef FESA_TESTS_THERMAL_GAUSS_ORACLE_H_", "#define FESA_TESTS_THERMAL_GAUSS_ORACLE_H_",
              "#include <array>", "namespace thermal_oracle {",
              "struct GaussPoint { double x_lo,x_hi,w_lo,w_hi; };"]
    for n in (12, 16):
        rows = rule(n)
        output.append(f"inline constexpr std::array<GaussPoint,{n}> kGauss{n}{{{{")
        output.extend("  {" + ",".join(x.hex() for x in row) + "}," for row in rows)
        output.append("}};")
    output.extend(["}", "#endif", ""])
    destination = Path(__file__).with_name("thermal_gauss_oracle.h")
    destination.write_text("\n".join(output), encoding="utf-8")
    print("Certified Legendre roots/weights: 12 + 16 distinct roots;", destination.name)
    fixtures = slab_fixtures()
    output = ["// Generated by thermal_oracle_generator.py. Outward binary64 brackets.",
              "#ifndef FESA_TESTS_THERMAL_SLAB_ORACLE_H_", "#define FESA_TESTS_THERMAL_SLAB_ORACLE_H_",
              "#include <array>", "namespace thermal_oracle {",
              "struct Bracket { double lo,hi; };",
              "struct SlabCase { double time; std::array<double,4> space; std::array<Bracket,13> fourier,continuous,be4096,be8192; };",
              "inline constexpr std::array<SlabCase,6> kSlab{{"]
    for time, space, *values in fixtures:
        output.append("{" + time.hex() + ",{{" + ",".join(x.hex() for x in space) + "}},")
        for rows in values:
            output.append("{{" + ",".join("{" + ",".join(x.hex() for x in bracket(row)) + "}" for row in rows) + "}},")
        output.append("},")
    output.extend(["}};", "}", "#endif", ""])
    Path(__file__).with_name("thermal_slab_oracle.h").write_text("\n".join(output), encoding="utf-8")