"""Separate exact polynomial cross-check; Python 3, standard library only.

Imports no audit_bs functions. Polynomials have integer coefficients; equality
compares every coefficient, not sampled floating-point evaluations. This verifies
algebraic identities, not empirical hypotheses or the missing original verifier.
Optional --output PATH saves the result, refusing to overwrite existing files.
"""
import argparse
import json
from pathlib import Path

NAMES = ('c', 's', 'dx', 'dy', 'dz', 'a', 'b', 'd', 'z')
ZERO = (0,) * len(NAMES)


class P:
    def __init__(self, terms=0):
        if isinstance(terms, int):
            terms = {ZERO: terms}
        self.terms = {m: a for m, a in terms.items() if a}

    @staticmethod
    def coerce(value):
        return value if isinstance(value, P) else P(value)

    def __add__(self, other):
        terms = dict(self.terms)
        for monomial, coefficient in self.coerce(other).terms.items():
            terms[monomial] = terms.get(monomial, 0) + coefficient
        return P(terms)

    __radd__ = __add__

    def __neg__(self):
        return P({m: -a for m, a in self.terms.items()})

    def __sub__(self, other):
        return self + (-self.coerce(other))

    def __rsub__(self, other):
        return self.coerce(other) + (-self)

    def __mul__(self, other):
        terms = {}
        for m, a in self.terms.items():
            for n, b in self.coerce(other).terms.items():
                power = tuple(x + y for x, y in zip(m, n))
                terms[power] = terms.get(power, 0) + a*b
        return P(terms)

    __rmul__ = __mul__

    def __eq__(self, other):
        return self.terms == self.coerce(other).terms


def variable(name):
    exponents = list(ZERO)
    exponents[NAMES.index(name)] = 1
    return P({tuple(exponents): 1})


def determinant(matrix):
    # Recursive cofactors, independently of the permutation implementation.
    if len(matrix) == 1:
        return P.coerce(matrix[0][0])
    return sum(((-1)**j)*matrix[0][j]*determinant(
        [row[:j] + row[j+1:] for row in matrix[1:]])
        for j in range(len(matrix)))


def calculate():
    c, s, dx, dy, dz, a, b, d, z = map(variable, NAMES)
    constraint = c*c + s*s - 1
    R = [[1,0,0,0], [0,1,0,0], [0,0,c,-s], [0,0,s,c]]
    checks = {}
    checks['rotation_orthogonal_mod_circle'] = all(
        sum(R[k][i]*R[k][j] for k in range(4)) - int(i == j)
        == (constraint if i == j and i >= 2 else P(0))
        for i in range(4) for j in range(4))
    checks['rotation_determinant_one_mod_circle'] = determinant(R)-1 == constraint
    checks['distance_identity_all_points_all_angles'] = (
        dx*dx + dy*dy + (c*dz)*(c*dz) + (s*dz)*(s*dz)
        - dx*dx - dy*dy - dz*dz == dz*dz*constraint)
    initial = [a, b, a+b+d, 0]
    checks['point_slice_hidden_coordinate'] = sum(R[3][i]*initial[i] for i in range(4)) == s*a+s*b+s*d
    V = determinant([[1,0,0], [0,1,0], [2,0,c]])
    W = determinant([[1,0,0], [0,1,0], [0,2,s]])
    checks['five_point_full_rank_circle_identity'] = V*V+W*W == c*c+s*s
    checks['point_tetra_full_dimension'] = determinant([[1,0,1], [0,1,1], [0,0,1]]) == 1
    endpoint_delta = [-z-z, P(0)-0]
    checks['diameter_is_visible_not_hidden'] = endpoint_delta == [-2*z, P(0)]
    return {'evidence_class': 'Exact polynomial identities; no physical measurements',
            'symbolic_checks': checks, 'all_passed': all(checks.values())}


if __name__ == '__main__':
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--output', type=Path)
    args = parser.parse_args()
    result = calculate()
    rendered = json.dumps(result, indent=2) + '\n'
    if args.output:
        with args.output.open('x', encoding='utf-8') as output:
            output.write(rendered)
    print(rendered, end='')
    raise SystemExit(0 if result['all_passed'] else 1)