#!/usr/bin/env python3
"""Exact checks for the supplementary field-interface note.

Standard library only. Explicit failures, so python3 -O does not disable them.
This is not a proof of the Galois correspondence, of the transcendence of pi,
or of a new cubature theorem.
"""

from __future__ import annotations

import json
import re
from fractions import Fraction
from math import gcd
from pathlib import Path

BITS = (
    (0, 0, 0),
    (1, 0, 0),
    (0, 1, 0),
    (0, 0, 1),
    (1, 1, 0),
    (1, 0, 1),
    (0, 1, 1),
    (1, 1, 1),
)
INDEX = {bits: i for i, bits in enumerate(BITS)}
QUADRATIC_RADICANDS = (2, 3, 5, 6, 10, 15, 30)
SCHEMA_ID = "omega-exact-coefficient/1"
PASSES = 0


def check(name: str, ok: bool) -> None:
    global PASSES
    if not ok:
        raise SystemExit(f"FAIL {name}")
    PASSES += 1
    print(f"PASS {name}")


def zero():
    return [Fraction(0) for _ in range(8)]


def mul(left, right):
    out = zero()
    for i, ai in enumerate(left):
        if ai == 0:
            continue
        e2i, e3i, e5i = BITS[i]
        for j, bj in enumerate(right):
            if bj == 0:
                continue
            e2, e3, e5 = e2i + BITS[j][0], e3i + BITS[j][1], e5i + BITS[j][2]
            scale = 1
            if e2 >= 2:
                scale *= 2
                e2 -= 2
            if e3 >= 2:
                scale *= 3
                e3 -= 2
            if e5 >= 2:
                scale *= 5
                e5 -= 2
            out[INDEX[(e2, e3, e5)]] += ai * bj * scale
    return out


def embed(element, signs):
    """Evaluate one embedding. signs are ±1 for sqrt(2), sqrt(3), sqrt(5)."""
    out = zero()
    for i, coeff in enumerate(element):
        e2, e3, e5 = BITS[i]
        sign = signs[0] ** e2 * signs[1] ** e3 * signs[2] ** e5
        out[i] = coeff * sign
    return out


def det(rows):
    size = len(rows)
    matrix = [list(row) for row in rows]
    sign = 1
    for col in range(size):
        pivot = next((row for row in range(col, size) if matrix[row][col] != 0), None)
        if pivot is None:
            return Fraction(0)
        if pivot != col:
            matrix[col], matrix[pivot] = matrix[pivot], matrix[col]
            sign = -sign
        pivot_value = matrix[col][col]
        for row in range(col + 1, size):
            if matrix[row][col] == 0:
                continue
            factor = matrix[row][col] / pivot_value
            for k in range(col, size):
                matrix[row][k] -= factor * matrix[col][k]
    value = Fraction(sign)
    for i in range(size):
        value *= matrix[i][i]
    return value


def element_in_f(u, v):
    """u + v*phi, phi=(1+sqrt(5))/2, as an element of K."""
    out = zero()
    out[0] = Fraction(u) + Fraction(v, 2)
    out[3] = Fraction(v, 2)
    return out


def golden_norm(u, v):
    u, v = Fraction(u), Fraction(v)
    return u * u + u * v - v * v


def field_norm(element):
    product = zero()
    product[0] = Fraction(1)
    for s2 in (1, -1):
        for s3 in (1, -1):
            for s5 in (1, -1):
                product = mul(product, embed(element, (s2, s3, s5)))
    return product


def odd_double_factorial(n):
    if n < -1:
        raise SystemExit("odd double factorial is defined here for n >= -1")
    if n <= 0:
        return 1
    value = 1
    for k in range(1, n + 1, 2):
        value *= k
    return value


def spherical_moment(i, j, k):
    """Normalized integral of x^{2i} y^{2j} z^{2k} on the unit sphere."""
    numerator = odd_double_factorial(2 * i - 1) * odd_double_factorial(2 * j - 1) * odd_double_factorial(2 * k - 1)
    denominator = odd_double_factorial(2 * i + 2 * j + 2 * k + 1)
    return Fraction(numerator, denominator)


def basis_label(bits):
    e2, e3, e5 = bits
    radicand = 1
    if e2:
        radicand *= 2
    if e3:
        radicand *= 3
    if e5:
        radicand *= 5
    if radicand == 1:
        return "1"
    return f"sqrt({radicand})"


def load_coefficient_schema():
    path = Path(__file__).with_name("coefficient_record.schema.json")
    if not path.is_file():
        raise SystemExit(f"FAIL missing coefficient schema {path.name}")
    try:
        schema = json.loads(path.read_text(encoding="utf-8"))
    except json.JSONDecodeError as exc:
        raise SystemExit(f"FAIL coefficient schema is not JSON: {exc}") from None
    if not isinstance(schema, dict):
        raise SystemExit("FAIL coefficient schema root is not an object")
    return schema


def resolve_schema(schema, root):
    if "$ref" not in schema:
        return schema
    ref = schema["$ref"]
    if not isinstance(ref, str) or not ref.startswith("#/"):
        raise SystemExit(f"FAIL unsupported schema ref {ref}")
    node = root
    for part in ref[2:].split("/"):
        part = part.replace("~1", "/").replace("~0", "~")
        if not isinstance(node, dict) or part not in node:
            raise SystemExit(f"FAIL broken schema ref {ref}")
        node = node[part]
    if not isinstance(node, dict):
        raise SystemExit(f"FAIL schema ref {ref} is not an object")
    return node


def schema_matches(instance, schema, root):
    """Evaluate the JSON Schema keywords this record uses. Not a general validator."""
    schema = resolve_schema(schema, root)
    if "allOf" in schema:
        if not all(schema_matches(instance, sub, root) for sub in schema["allOf"]):
            return False
    expected = schema.get("type")
    if expected == "object":
        if not isinstance(instance, dict) or isinstance(instance, list):
            return False
        props = schema.get("properties", {})
        for key in schema.get("required", []):
            if key not in instance:
                return False
        if schema.get("additionalProperties") is False and any(key not in props for key in instance):
            return False
        for key, sub in props.items():
            if key in instance and not schema_matches(instance[key], sub, root):
                return False
    elif expected == "array":
        if not isinstance(instance, list):
            return False
        if "minItems" in schema and len(instance) < schema["minItems"]:
            return False
        if "maxItems" in schema and len(instance) > schema["maxItems"]:
            return False
        item_schema = schema.get("items")
        if isinstance(item_schema, dict) and any(not schema_matches(item, item_schema, root) for item in instance):
            return False
    elif expected == "string":
        if not isinstance(instance, str):
            return False
    elif expected is not None:
        raise SystemExit(f"FAIL schema uses an unsupported type {expected}")
    if "const" in schema and instance != schema["const"]:
        return False
    if "pattern" in schema:
        if not isinstance(instance, str) or "\n" in instance or "\r" in instance:
            return False
        if re.search(schema["pattern"], instance) is None:
            return False
    return True


def schema_document_enforces_format(schema):
    if schema.get("$schema") != "https://json-schema.org/draft/2020-12/schema":
        return False
    if schema.get("type") != "object" or schema.get("additionalProperties") is not False:
        return False
    if schema.get("required") != ["schema", "coefficients"]:
        return False
    props = schema.get("properties")
    if not isinstance(props, dict):
        return False
    if props.get("schema", {}).get("const") != SCHEMA_ID:
        return False
    coefficients = props.get("coefficients")
    if not isinstance(coefficients, dict):
        return False
    if coefficients.get("type") != "array" or coefficients.get("minItems") != 8 or coefficients.get("maxItems") != 8:
        return False
    if coefficients.get("items", {}).get("$ref") != "#/$defs/coefficient":
        return False
    coefficient = schema.get("$defs", {}).get("coefficient")
    if not isinstance(coefficient, dict):
        return False
    if coefficient.get("type") != "object" or coefficient.get("additionalProperties") is not False:
        return False
    if coefficient.get("required") != ["n", "d"]:
        return False
    fields = coefficient.get("properties")
    if not isinstance(fields, dict):
        return False
    if fields.get("n", {}).get("type") != "string" or fields.get("d", {}).get("type") != "string":
        return False
    if fields.get("n", {}).get("pattern") != r"^(0|-?[1-9][0-9]*)$":
        return False
    if fields.get("d", {}).get("pattern") != r"^[1-9][0-9]*$":
        return False
    return True


def coefficient_record(coefficients, schema_id=SCHEMA_ID):
    return {"schema": schema_id, "coefficients": coefficients}


def eight_coefficients(replace=None, index=0):
    coefficients = [{"n": "0", "d": "1"} for _ in range(8)]
    if replace is not None:
        coefficients[index] = replace
    return coefficients


def is_reduced_coefficient(coefficient):
    if not isinstance(coefficient, dict):
        return False
    numerator = coefficient.get("n")
    denominator = coefficient.get("d")
    if not isinstance(numerator, str) or not isinstance(denominator, str):
        return False
    if re.fullmatch(r"(0|-?[1-9][0-9]*)", numerator) is None:
        return False
    if re.fullmatch(r"[1-9][0-9]*", denominator) is None:
        return False
    n = int(numerator)
    d = int(denominator)
    if gcd(abs(n), d) != 1:
        return False
    if n == 0 and (numerator != "0" or denominator != "1"):
        return False
    return True


def reduced_record(record):
    if not isinstance(record, dict):
        return False
    coefficients = record.get("coefficients")
    if not isinstance(coefficients, list) or len(coefficients) != 8:
        return False
    return all(is_reduced_coefficient(item) for item in coefficients)


def main():
    sign_rows = []
    for s2 in (1, -1):
        for s3 in (1, -1):
            for s5 in (1, -1):
                sign_rows.append(
                    [
                        Fraction(s2 ** e2 * s3 ** e3 * s5 ** e5)
                        for e2, e3, e5 in BITS
                    ]
                )
    sign_det = det(sign_rows)
    check("eight basis embeddings have nonzero sign determinant", sign_det != 0)
    check("basis has eight ordered square-free monomials", len(BITS) == 8)
    check("F sits inside K because sqrt(5) is a basis element", INDEX[(0, 0, 1)] == 3)

    phi = element_in_f(0, 1)
    one = element_in_f(1, 0)
    phi_square = element_in_f(1, 1)
    check("phi = (1+sqrt(5))/2 uses only the rational and sqrt(5) slots", phi == [Fraction(1, 2), 0, 0, Fraction(1, 2), 0, 0, 0, 0])
    check("1 and phi squared are distinct and share the golden norm 1", one != phi_square and golden_norm(1, 0) == 1 and golden_norm(1, 1) == 1)
    check("1 and phi share the larger norm and not the golden norm", golden_norm(1, 0) == 1 and golden_norm(0, 1) == -1)

    samples = ((0, 1), (1, 0), (1, 1), (2, -3), (Fraction(3, 2), Fraction(-4, 5)))
    for u, v in samples:
        got = field_norm(element_in_f(u, v))
        expected = zero()
        expected[0] = golden_norm(u, v) ** 4
        check(f"larger norm is the fourth power of the golden norm at ({u}, {v})", got == expected)

    check("13 is not one of the seven square-free radicands of the quadratic subfields", 13 not in QUADRATIC_RADICANDS)

    n = (Fraction(2), Fraction(3), Fraction(0))
    gram = n[0] * n[0] + n[1] * n[1] + n[2] * n[2]
    check("n=(2,3,0) has squared length 13", gram == 13)
    reflection = [
        [Fraction(5, 13), Fraction(-12, 13), Fraction(0)],
        [Fraction(-12, 13), Fraction(-5, 13), Fraction(0)],
        [Fraction(0), Fraction(0), Fraction(1)],
    ]
    householder = [[Fraction(1 if i == j else 0) for j in range(3)] for i in range(3)]
    for i in range(3):
        for j in range(3):
            householder[i][j] -= 2 * n[i] * n[j] / gram
    check("displayed reflection equals I-2nn^T/(n^T n)", reflection == householder)

    def matmul(left, right):
        return [
            [sum(left[i][k] * right[k][j] for k in range(3)) for j in range(3)]
            for i in range(3)
        ]

    def apply(matrix, vector):
        return [sum(matrix[i][j] * vector[j] for j in range(3)) for i in range(3)]

    identity = [[Fraction(1 if i == j else 0) for j in range(3)] for i in range(3)]
    transpose = [[reflection[j][i] for j in range(3)] for i in range(3)]
    check("R squared is the identity", matmul(reflection, reflection) == identity)
    check("R is orthogonal", matmul(transpose, reflection) == identity)
    check("Rn = -n", apply(reflection, list(n)) == [-n[0], -n[1], -n[2]])

    def disagree_at(left, right):
        hits = [d for d, (a, b) in enumerate(zip(left, right)) if a != b]
        return min(hits) if hits else None

    series_a = (0, 1, 2, 3)
    series_b = (0, 1, 2, 4)
    star = disagree_at(series_a, series_b)
    check(
        "integer first disagreement is attained",
        star == 3 and series_a[:star] == series_b[:star] and series_a[star] != series_b[star],
    )

    def continuous_pair(t):
        return Fraction(0), (Fraction(0) if t <= 0 else t)

    infimum = Fraction(0)
    left0, right0 = continuous_pair(infimum)
    left_pos, right_pos = continuous_pair(Fraction(1, 10))
    check(
        "a continuous infimum need not be a disagreement point",
        left0 == right0 and left_pos != right_pos,
    )

    y, scale = 10, 3
    quotient, remainder = divmod(y, scale)
    check("positive floor identity recovers the input", scale * quotient + remainder == y and 0 <= remainder < scale)

    moment = spherical_moment(2, 0, 1) - spherical_moment(1, 1, 1) + spherical_moment(0, 2, 1)
    check("normalized integral of J squared is 1/21", moment == Fraction(1, 21) and moment > 0)
    check("x^2 y^2 z^2 has normalized spherical mean 1/105", spherical_moment(1, 1, 1) == Fraction(1, 105))

    schema = load_coefficient_schema()
    check(
        "coefficient order is the displayed basis order recorded in basis_order and BITS",
        schema.get("basis_order") == [basis_label(bits) for bits in BITS],
    )
    check("coefficient schema enforces the omega-exact-coefficient/1 record format", schema_document_enforces_format(schema))

    def accepted(record):
        return schema_matches(record, schema, schema) and reduced_record(record)

    phi_coefficients = eight_coefficients()
    phi_coefficients[0] = {"n": "1", "d": "2"}
    phi_coefficients[3] = {"n": "1", "d": "2"}
    check("a reduced eight-coefficient record is accepted", accepted(coefficient_record(phi_coefficients)))

    missing_records = [
        {"schema": SCHEMA_ID},
        coefficient_record(eight_coefficients()[:7]),
        coefficient_record(eight_coefficients() + [{"n": "0", "d": "1"}]),
        coefficient_record(eight_coefficients({"d": "1"})),
        coefficient_record(eight_coefficients({"n": "1"})),
        coefficient_record(eight_coefficients("0/1")),
        {"schema": SCHEMA_ID, "coefficients": eight_coefficients(), "trace": "0"},
        coefficient_record(eight_coefficients({"n": "0", "d": "1", "norm": "1"})),
        coefficient_record(eight_coefficients(), schema_id="omega-exact-coefficient/2"),
    ]
    check(
        "schema rejects missing coefficients and records that are not eight coefficient objects",
        all(not schema_matches(record, schema, schema) for record in missing_records),
    )

    malformed = [
        {"n": "01", "d": "1"},
        {"n": "-0", "d": "1"},
        {"n": "+1", "d": "1"},
        {"n": "1.5", "d": "1"},
        {"n": "", "d": "1"},
        {"n": "1 ", "d": "1"},
        {"n": " 1", "d": "1"},
        {"n": "1a", "d": "1"},
        {"n": "1\n", "d": "1"},
        {"n": "1", "d": "01"},
        {"n": "1", "d": "-1"},
        {"n": "1", "d": "+2"},
        {"n": "1", "d": "1.0"},
        {"n": 1, "d": "1"},
        {"n": "1", "d": 1},
    ]
    check(
        "schema rejects malformed integer strings",
        all(not schema_matches(coefficient_record(eight_coefficients(item)), schema, schema) for item in malformed),
    )

    zero_denominators = [
        {"n": "1", "d": "0"},
        {"n": "0", "d": "0"},
        {"n": "1", "d": "00"},
        {"n": "1", "d": "-0"},
    ]
    check(
        "schema rejects a zero denominator",
        all(not schema_matches(coefficient_record(eight_coefficients(item)), schema, schema) for item in zero_denominators),
    )

    unreduced = [
        {"n": "2", "d": "4"},
        {"n": "-2", "d": "4"},
        {"n": "3", "d": "6"},
        {"n": "-4", "d": "6"},
        {"n": "0", "d": "2"},
        {"n": "0", "d": "10"},
    ]
    reduced_samples = [
        {"n": "0", "d": "1"},
        {"n": "-1", "d": "2"},
        {"n": "2", "d": "1"},
        {"n": "-4", "d": "5"},
    ]
    check(
        "reduced-form validator rejects unreduced fractions and any zero other than 0/1",
        all(schema_matches(coefficient_record(eight_coefficients(item)), schema, schema) for item in unreduced)
        and all(not is_reduced_coefficient(item) for item in unreduced)
        and all(is_reduced_coefficient(item) for item in reduced_samples)
        and all(not accepted(coefficient_record(eight_coefficients(item))) for item in unreduced),
    )

    print(f"TOTAL {PASSES} PASS")
    print("Checked: basis sign determinant, golden-to-larger norm, reflection, disagreement, floor identity, J^2 moment, coefficient schema.")
    print("Not checked here: transcendence of pi, a general Galois-correspondence proof, or a new cubature rule.")


if __name__ == "__main__":
    main()