"""Exact ten-sample reconstruction on H_{<=2} + R*J.

Python 3.10+; standard library only. Run normally or with python -O.
All comparisons use exact algebraic arithmetic, never a floating tolerance.

Field representation:
  Q(sqrt(5), r, sqrt(3)), r^2=(5+sqrt(5))/2=1+phi^2.
Basis bit order: sqrt(5), r, sqrt(3); coefficients are rational.
"""
from dataclasses import dataclass
from fractions import Fraction as F
from itertools import product
from math import sqrt


def require(ok, label):
    if not ok:
        raise RuntimeError(label)


@dataclass(frozen=True)
class Exact:
    coeffs: tuple

    def __post_init__(self):
        require(len(self.coeffs) == 8, 'eight algebraic coefficients')
        object.__setattr__(self, 'coeffs', tuple(map(F, self.coeffs)))

    @classmethod
    def rational(cls, value):
        return cls((F(value),) + (F(0),)*7)

    def __add__(self, other):
        other = coerce(other)
        return Exact(tuple(a+b for a,b in zip(self.coeffs, other.coeffs)))

    __radd__ = __add__

    def __neg__(self):
        return Exact(tuple(-a for a in self.coeffs))

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

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

    def __mul__(self, other):
        other = coerce(other)
        out = [F(0)]*8
        for i,a in enumerate(self.coeffs):
            if not a:
                continue
            for j,b in enumerate(other.coeffs):
                if not b:
                    continue
                e5 = (i & 1) + (j & 1)
                er = ((i >> 1) & 1) + ((j >> 1) & 1)
                e3 = ((i >> 2) & 1) + ((j >> 2) & 1)
                factor = a*b * (3 if e3 == 2 else 1)
                e3 %= 2
                # r^2 = 5/2 + sqrt(5)/2.
                terms = [(e5,F(5,2)), (e5+1,F(1,2))] if er == 2 else [(e5,F(1))]
                er %= 2
                for ep,cp in terms:
                    v = factor*cp * 5**(ep//2)
                    index = ep%2 + 2*er + 4*e3
                    out[index] += v
        return Exact(tuple(out))

    __rmul__ = __mul__

    def __truediv__(self, value):
        # General inverses are not needed for this explicit decoder.
        require(isinstance(value,(int,F)) and value != 0, 'rational divisor')
        return Exact(tuple(c/F(value) for c in self.coeffs))

    def __pow__(self, n):
        require(isinstance(n,int) and n >= 0, 'nonnegative integer exponent')
        out = coerce(1)
        for _ in range(n):
            out = out*self
        return out

    def __bool__(self):
        return any(self.coeffs)

    def decimal(self):
        rr = sqrt((5+sqrt(5))/2)
        return sum(float(c) * sqrt(5)**(i&1) * rr**((i>>1)&1)
                   * sqrt(3)**((i>>2)&1) for i,c in enumerate(self.coeffs))


def coerce(value):
    return value if isinstance(value,Exact) else Exact.rational(value)


def radical(index):
    c = [F(0)]*8
    c[index] = F(1)
    return Exact(tuple(c))


ZERO, ONE = coerce(0), coerce(1)
ROOT5, RADIUS, ROOT3 = radical(1), radical(2), radical(4)
PHI = (ONE+ROOT5)/2
PHI_PRIME = ONE-PHI
R_INV = RADIUS*(5-ROOT5)/10
ROOT3_INV = ROOT3/3
ROOT5_INV = ROOT5/5


def dot(a,b):
    return sum((x*y for x,y in zip(a,b)),ZERO)


def vscale(c,v):
    return tuple(c*x for x in v)


def matvec(m,v):
    return tuple(dot(row,v) for row in m)


def transpose(m):
    return tuple(zip(*m))


def matmul(a,b):
    return tuple(tuple(dot(row,col) for col in transpose(b)) for row in a)


AXES_RAW = [(0,1,PHI),(1,PHI,0),(PHI,0,1),
            (0,1,-PHI),(1,-PHI,0),(PHI,0,-1)]
AXES = [tuple(coerce(x)*R_INV for x in v) for v in AXES_RAW]
SAMPLE_NODES = [v for a in AXES[:3] for v in (a,vscale(-1,a))] + AXES[3:]
PSTAR = (ROOT3_INV,)*3
# inverse of rows a_1, a_2, a_3; 1/phi^2 = 2-phi.
B_INV_FACTOR = RADIUS*(2-PHI)/2
B_INV = tuple(tuple(B_INV_FACTOR*x for x in row) for row in
              [(-PHI,ONE,PHI**2),(ONE,PHI**2,-PHI),(PHI**2,-PHI,ONE)])

# xi=R*x are aligned coordinates. r_-^{-1}=r/sqrt(5).
ALIGN = ((R_INV,PHI*R_INV,ZERO),
         (RADIUS*ROOT5_INV,PHI_PRIME*RADIUS*ROOT5_INV,ZERO),
         (ZERO,ZERO,ONE))


# Sparse polynomials in aligned coordinates (X,Y,Z).
CONST = {(0,0,0):ONE}


def clean(p):
    return {e:coerce(c) for e,c in p.items() if coerce(c)}


def add(*polys):
    out = {}
    for p in polys:
        for e,c in p.items():
            out[e] = out.get(e,ZERO)+c
    return clean(out)


def scale(c,p):
    return clean({e:coerce(c)*v for e,v in p.items()})


def mul(p,q):
    out = {}
    for e,c in p.items():
        for f,d in q.items():
            g = tuple(a+b for a,b in zip(e,f))
            out[g] = out.get(g,ZERO)+c*d
    return clean(out)


def variable(i):
    e = [0]*3
    e[i] = 1
    return {tuple(e):ONE}


VARS = tuple(variable(i) for i in range(3))
R2 = add(*(mul(x,x) for x in VARS))
J = {(1,1,1):ROOT5}


def linear(v):
    return add(*(scale(c,x) for c,x in zip(v,VARS)))


def evaluate(p,node):
    return sum((c*node[0]**e[0]*node[1]**e[1]*node[2]**e[2]
                for e,c in p.items()),ZERO)


def derivative(p,i):
    out = {}
    for e,c in p.items():
        if e[i]:
            f = list(e)
            f[i] -= 1
            out[tuple(f)] = c*e[i]
    return clean(out)


def laplacian(p):
    return add(*(derivative(derivative(p,i),i) for i in range(3)))


EVEN = [add(scale(F(1,6),CONST),
            scale(F(5,4),add(mul(linear(a),linear(a)),scale(F(-1,3),R2))))
        for a in AXES]
J_STAR = ROOT5*ROOT3/9
J_STAR_INV = 3*ROOT3*ROOT5/5


def decode(values):
    require(len(values) == 10, 'exactly ten ordered samples')
    values = tuple(map(coerce,values))
    differences = tuple((values[2*j]-values[2*j+1])/2 for j in range(3))
    ell = matvec(B_INV,differences)
    masses = [(values[2*j]+values[2*j+1])/2 for j in range(3)]
    masses += [values[6+j]-dot(ell,AXES[3+j]) for j in range(3)]
    u = add(linear(ell),*(scale(c,e) for c,e in zip(masses,EVEN)))
    beta = (values[9]-evaluate(u,PSTAR))*J_STAR_INV
    return add(u,scale(beta,J)), u, beta


def compose_linear(poly, matrix):
    substitution = [linear(row) for row in matrix]
    out = {}
    for e,c in poly.items():
        term = scale(c,CONST)
        for i,power in enumerate(e):
            for _ in range(power):
                term = mul(term,substitution[i])
        out = add(out,term)
    return out


def physical_j(point):
    x,y,z = point
    return z*(x*x+x*y-y*y)


def exact_text(value):
    labels = ['1','sqrt5','r','sqrt5*r','sqrt3','sqrt5*sqrt3','r*sqrt3','sqrt5*r*sqrt3']
    terms = []
    for c,label in zip(value.coeffs,labels):
        if c:
            terms.append(str(c) if label=='1' else '('+str(c)+')*'+label)
    return ' + '.join(terms) if terms else '0'


def run():
    passed = []
    identity = tuple(tuple(ONE if i==j else ZERO for j in range(3)) for i in range(3))
    require(ROOT5**2==coerce(5) and ROOT3**2==coerce(3), 'radical relations')
    require(RADIUS**2==ONE+PHI**2 and RADIUS*R_INV==ONE, 'radius and inverse')
    require(J_STAR*J_STAR_INV==ONE, 'extra sample denominator')
    passed.append('algebraic constants and nonzero extra-sample denominator')
    require(matmul(AXES[:3],B_INV)==identity, 'explicit inverse')
    require(matmul(B_INV,AXES[:3])==identity, 'inverse in both directions')
    passed.append('explicit 3 by 3 inverse, both directions')
    require(all(dot(a,a)==ONE for a in AXES), 'unit axes')
    require(all(dot(AXES[i],AXES[j])**2==coerce(F(1,5))
                for i in range(6) for j in range(i)), 'equiangular axes')
    passed.append('six unit axes and all 15 squared cross-inner-products')
    require(all(evaluate(EVEN[j],AXES[i])==(ONE if i==j else ZERO)
                for i in range(6) for j in range(6)), 'even cardinality')
    require(all(not laplacian(e) for e in EVEN), 'harmonic even basis')
    passed.append('36 even-cardinal identities and six harmonicity identities')
    require(all(evaluate(e,PSTAR)==PHI/3 for e in EVEN[:3]), 'positive star values')
    require(all(evaluate(e,PSTAR)==PHI_PRIME/3 for e in EVEN[3:]), 'negative star values')
    passed.append('six simplified even-interpolant values at the extra point')
    require(all(evaluate(J,n)==ZERO for n in SAMPLE_NODES), 'invisible cubic')
    require(evaluate(J,PSTAR)==J_STAR and not laplacian(J), 'harmonic cubic and star')
    passed.append('nine cubic zeros, harmonicity, and exact extra sample')
    cardinals = []
    for k in range(10):
        v = [ONE if i==k else ZERO for i in range(10)]
        cardinals.append(decode(v)[0])
    nodes = SAMPLE_NODES+[PSTAR]
    require(all(evaluate(cardinals[j],nodes[i])==(ONE if i==j else ZERO)
                for i in range(10) for j in range(10)), 'ten-point cardinal matrix')
    passed.append('100 exact ten-point cardinal identities')
    # Independently assemble the closed cardinal formula stated in the proof.
    wpoly = [add(VARS[i],*(scale(-AXES[k][i],EVEN[k]) for k in range(3,6)))
             for i in range(3)]
    tpoly = [add(*(scale(B_INV[i][j],wpoly[i]) for i in range(3))) for j in range(3)]
    low_cardinals = [scale(F(1,2),add(EVEN[j],scale(sign,tpoly[j])))
                     for j in range(3) for sign in (1,-1)] + EVEN[3:]
    cubic_cardinal = scale(J_STAR_INV,J)
    closed_cardinals = [add(l,scale(-evaluate(l,PSTAR),cubic_cardinal))
                        for l in low_cardinals] + [cubic_cardinal]
    require(closed_cardinals==cardinals, 'closed cardinal formula equals staged decoder')
    passed.append('ten closed cardinal formulas equal the staged decoder')
    basis = [CONST,*VARS,mul(VARS[0],VARS[1]),mul(VARS[0],VARS[2]),
             mul(VARS[1],VARS[2]),add(mul(VARS[0],VARS[0]),scale(-1,mul(VARS[1],VARS[1]))),
             add(mul(VARS[1],VARS[1]),scale(-1,mul(VARS[2],VARS[2]))),J]
    for p in basis:
        recovered,_,_ = decode([evaluate(p,n) for n in nodes])
        require(recovered==p, 'polynomial reconstruction identity')
    passed.append('ten exact polynomial reconstruction identities')
    require(matmul(ALIGN,transpose(ALIGN))==identity, 'orthogonal alignment')
    require(compose_linear(J,ALIGN)=={(2,0,1):ONE,(1,1,1):ONE,(0,2,1):-ONE},
            'J equals sqrt(5)XYZ under alignment')
    passed.append('orthogonal alignment and physical cubic identity')
    physical_nodes = [matvec(transpose(ALIGN),n) for n in nodes]
    def physical_example(n):
        x,y,z = n
        return 2+3*x+x*x-y*y+4*physical_j(n)
    result,_,beta = decode([physical_example(n) for n in physical_nodes])
    probe = tuple(map(coerce,[F(1,5),F(-3,10),F(1,3)]))
    result_at_probe = evaluate(result,matvec(ALIGN,probe))
    require(beta==coerce(4) and result_at_probe==coerce(F(721,300)), 'worked example')
    passed.append('physical worked example: beta=4 and probe value=721/300')
    physical_poly = compose_linear(result,ALIGN)
    # These are exact finite transport checks, supplementing the general proof.
    for tau in [F(1,2),F(2,3)]:
        A = ((0,tau,0),(tau,tau,0),(0,0,-1/tau**2))
        Ainv = ((-1/tau,1/tau,0),(1/tau,0,0),(0,0,-tau**2))
        A = tuple(tuple(map(coerce,row)) for row in A)
        Ainv = tuple(tuple(map(coerce,row)) for row in Ainv)
        transported_poly = compose_linear(physical_poly,Ainv)
        measurements = [evaluate(transported_poly,matvec(A,n)) for n in physical_nodes]
        transported_result,_,transported_beta = decode(measurements)
        require(transported_result==result and transported_beta==coerce(4), 'transported values')
    passed.append('exact sample transport at two rational scales')
    try:
        decode([0]*9)
    except RuntimeError:
        pass
    else:
        raise RuntimeError('missing tenth sample was accepted')
    passed.append('wrong sample count rejected')
    for item in passed:
        print('PASS',item)
    print('TOTAL',len(passed),'PASS')
    print('EXAMPLE beta = 4; f(1/5,-3/10,1/3) = 721/300')
    print('All equality checks above are exact. No floating-point acceptance tests.')
    print('EXACT INPUTS: sqrt5=sqrt(5), sqrt3=sqrt(3), r=sqrt((5+sqrt(5))/2)')
    labels = ['v1+','v1-','v2+','v2-','v3+','v3-','v4','v5','v6','v*']
    for label,n in zip(labels,physical_nodes):
        print(label,'=',exact_text(physical_example(n)))


if __name__=='__main__':
    run()