"""OMEGA transition invariants: exact algebraic witnesses and finite checks.

Run: python3 -B OMEGA_Transition_Invariants_Proof.py
Run: python3 -B -O OMEGA_Transition_Invariants_Proof.py
Dependencies: Python standard library only, Python 3.10 or later.

The polynomial engine normalizes expressions over Q[phi]/(phi^2-phi-1).
Other symbols are independent commuting variables. pi remains an exact symbol.
Zero normal forms establish the displayed polynomial identities for all
substitutions satisfying the declared assumptions. Rational identities use
cross multiplication and require every displayed denominator to be nonzero.

ASSUMPTIONS
* phi=(1+sqrt(5))/2. Geometry uses fixed s>0, h>0.
* G=h*phi and C=h*pi name the declared selected paths, not the whole outline.
* Similarities preserve the shape. Projective maps have nonzero determinant,
  no pole at a selected mark, and four distinct labeled marks.
* A resolution change keeps the same polynomial field and recalibrates its
  coefficients. Raw readings at different physical locations need not agree.
* The wave proof uses u_tt=u_xx, speed one, no forcing, and degree<=5 initial
  field and velocity. The moving paths are x=x0-t and x=x0+t.
* Q,S are positive integers, Q<S, and the ledger retains every ordered residual.
* RN scaling compares the SAME solution. h>0, scale factor c>0, m!=0.
  The normalized Einstein residual uses natural units hbar=c_light=1 and
  orthonormal components at one fixed physical point.

These are algebraic witnesses plus implementation tests. They are not a Lean
certificate, an independent review, a physical experiment, or a unification
proof. Exact algebra does not establish noise stability or real calibration.
The canonical sine is not evaluated exactly here. Geometry identities retain
symbolic s; finite sample failures use s=1/2, a declared family member.
"""

from fractions import Fraction as F
from math import factorial
import sys


def fib(n):
    a,b=0,1
    for _ in range(n):
        a,b=b,a+b
    return a


class P:
    """Sparse multivariate polynomial with exact rational coefficients."""
    def __init__(self, terms=None):
        self.terms={}
        for monomial,coefficient in (terms or {}).items():
            coefficient=F(coefficient)
            powers=dict(monomial)
            n=powers.pop('phi',0)
            base=tuple(sorted((k,v) for k,v in powers.items() if v))
            # phi^n = F_(n-1) + F_n*phi for n>=1.
            expansions=[(0,1)] if n==0 else [(0,fib(n-1)),(1,fib(n))]
            for degree,multiplier in expansions:
                key=tuple(sorted(base+(('phi',degree),))) if degree else base
                self.terms[key]=self.terms.get(key,F(0))+coefficient*multiplier
        self.terms={k:v for k,v in self.terms.items() if v}

    @staticmethod
    def coerce(x):
        if isinstance(x,P): return x
        if type(x) not in (int,F): raise TypeError('use exact integers or fractions')
        return P({():F(x)})

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

    __radd__=__add__

    def __neg__(self): return P({k:-v for k,v 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 ka,va in self.terms.items():
            for kb,vb in self.coerce(other).terms.items():
                powers=dict(ka)
                for name,power in kb: powers[name]=powers.get(name,0)+power
                key=tuple(sorted(powers.items()))
                terms[key]=terms.get(key,F(0))+va*vb
        return P(terms)

    __rmul__=__mul__

    def __pow__(self,n):
        if type(n) is not int or n<0: raise ValueError('nonnegative integer power required')
        result=P.coerce(1)
        for _ in range(n): result*=self
        return result

    def derivative(self,name,order=1):
        if type(order) is not int or order<0: raise ValueError('invalid derivative order')
        result=self
        for _ in range(order):
            terms={}
            for key,value in result.terms.items():
                powers=dict(key); degree=powers.get(name,0)
                if degree:
                    powers[name]=degree-1
                    reduced=tuple(sorted((k,v) for k,v in powers.items() if v))
                    terms[reduced]=terms.get(reduced,F(0))+value*degree
            result=P(terms)
        return result

    def substitute(self,**replacements):
        result=P.coerce(0)
        for key,value in self.terms.items():
            term=P.coerce(value)
            for name,power in key:
                term*=P.coerce(replacements[name])**power if name in replacements else var(name)**power
            result+=term
        return result


def var(name): return P({((name,1),):F(1)})


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


def identity(expression,label):
    require(not expression.terms,label+' failed')


def different(expression,label):
    require(bool(expression.terms),label+' failed to expose a difference')


def ratio_identity(numerator,denominator,reference_numerator,reference_denominator,label):
    identity(numerator*reference_denominator-reference_numerator*denominator,label)


def evolve_polynomial(p,v,x,t):
    """Exact finite d'Alembert series for any supplied polynomial initial data."""
    degree=max([dict(k).get('x',0) for k in p.terms|v.terms] or [0])
    u=P.coerce(0)
    for n in range(degree//2+1):
        u+=p.derivative('x',2*n)*t**(2*n)*F(1,factorial(2*n))
        u+=v.derivative('x',2*n)*t**(2*n+1)*F(1,factorial(2*n+1))
    return u


def main():
    print('OMEGA transition invariant witnesses v1.0')
    print('Optimization mode:',sys.flags.optimize)
    phi,s,h,pi,c=map(var,['phi','s','h','pi','scale'])
    x,t,x0=map(var,['x','t','x0'])

    # Engine checks distinguish formal variables from the golden extension.
    identity(phi**2-phi-1,'golden identity')
    identity(phi**7-(13*phi+8),'higher golden power')
    identity(((x+2)*(x-3)).derivative('x')-(2*x-1),'polynomial derivative')
    identity((x*x+t).substitute(x=x0-t)-((x0-t)**2+t),'polynomial substitution')
    different(x-1,'nonzero polynomial detection')
    print('PASS 1: exact polynomial-engine consistency checks')

    G,C=h*phi,h*pi
    identity(G+C-h*(phi+pi),'Omega0 calibration')
    ratio_identity(G,G+C,phi,phi+pi,'lambda selected-path ratio')
    A,B=s+phi**2,s+1
    z=(-A,-B,-s,s,B,A)
    expected=(phi,P.coerce(1),2*s,P.coerce(1),phi)
    origin=var('origin')
    marks=tuple(origin+h*zi for zi in z)
    for j,gap in enumerate(expected): identity(marks[j+1]-marks[j]-h*gap,'normalized gap')
    co,si,dx,dy=map(var,['cosine','sine','dx','dy'])
    identity((co*dx-si*dy)**2+(si*dx+co*dy)**2
             -(co**2+si**2)*(dx**2+dy**2),'rotation norm factor')
    print('PASS 2: calibration, five normalized gaps, and rotation norm identity')

    aa,bb,cc,dd=map(var,['a','b','c','d'])
    def diff_numerator(u,v): return (aa*u+bb)*(cc*v+dd)-(aa*v+bb)*(cc*u+dd)
    identity(diff_numerator(x,x0)-(aa*dd-bb*cc)*(x-x0),'projective difference factor')
    x1,x2,x3,x4=map(var,['x1','x2','x3','x4'])
    ratio_identity(diff_numerator(x3,x1)*diff_numerator(x4,x2),
                   diff_numerator(x3,x2)*diff_numerator(x4,x1),
                   (x3-x1)*(x4-x2),(x3-x2)*(x4-x1),'projective cross-ratio')
    identity((A+B)**2-(4*A*B+phi**2),'OMEGA reference cross-ratio')
    print('PASS 3: general projective cross-ratio and OMEGA specialization')

    hp=var('new_h')
    for j in range(6):
        coefficient=var('coefficient_'+str(j))
        ratio_identity(coefficient*h**j,h**j,coefficient*hp**j,hp**j,
                       'same physical coefficient at both resolutions')
    print('PASS 4: exact recovery of all six coefficients after reference-length changes')

    p=sum(var('p'+str(j))*x**j for j in range(6))
    v=sum(var('v'+str(j))*x**j for j in range(6))
    u=evolve_polynomial(p,v,x,t)
    identity(u.derivative('t',2)-u.derivative('x',2),'wave equation')
    identity(u.substitute(t=0)-p,'initial field')
    identity(u.derivative('t').substitute(t=0)-v,'initial velocity')
    plus=u.derivative('t')+u.derivative('x')
    minus=u.derivative('t')-u.derivative('x')
    identity(plus.substitute(x=x0-t)-(v+p.derivative('x')).substitute(x=x0),
             'left-moving characteristic invariant')
    identity(minus.substitute(x=x0+t)-(v-p.derivative('x')).substitute(x=x0),
             'right-moving characteristic invariant')
    print('PASS 5: both characteristic invariants for arbitrary P5 coefficients, x0, and time')

    q,ss,d,y,numerator,weight=map(var,['Q','S','state','next_state','prefix_numerator','S_to_N'])
    rho=q*d-ss*y
    identity(ss*y+rho-q*d,'one-step reconstruction')
    identity(d+0-d,'ledger base case N=0')
    # Prefix recurrence W_(N+1)=Q W_N+rho_N S^N, B_N=S^N d_N+W_N.
    next_numerator=q*numerator+weight*(ss*y+rho-q*d)
    identity(next_numerator-q*numerator,'ledger induction step')
    Q,S=339949771344778,10**15
    comparisons=0;max_hit=0
    for seed in range(32768):
        state=seed; W=0; qpower=spower=1; hit=0 if seed==0 else None
        for n in range(1,24):
            next_state,residual=divmod(Q*state,S)
            W=Q*W+residual*spower
            qpower*=Q;spower*=S
            require(spower*next_state+W==qpower*seed,'finite prefix invariant')
            require((S*next_state+residual)//Q==state,'inverse step')
            state=next_state
            if hit is None and state==0: hit=n
            comparisons+=1
        require(hit is not None,'finite extinction domain')
        max_hit=max(max_hit,hit)
    require(max_hit==10,'maximum canonical extinction time')
    print('PASS 6: ledger base/induction algebra and',comparisons,'canonical prefix invariants')

    m,chi=map(var,['mass_length','charge_squared_length'])
    h2=h**2; newh2=(c*h)**2
    # (chi/h^2)/(m^2/h^2), retaining its nonzero denominators.
    ratio_identity(chi*h2,h2*m**2,chi*newh2,newh2*m**2,'RN ruler-change ratio')
    ratio_identity(chi*h2,h2*m**2,chi,m**2,'RN physical ratio')
    g,stress,kappa,lam,eta=map(var,['Gcomponent','Tcomponent','kappa','Lambda','metric_component'])
    E=g+lam*eta-kappa*stress
    identity((c*h)**2*E-c**2*(h**2*E),'normalized Einstein residual scaling')
    print('PASS 7: RN dimensionless invariant and normalized residual scaling')

    # Negative controls: the code must not label these quantities invariant.
    different((c*G-G).substitute(scale=2,h=1),'raw length changes under rescaling')
    old_s=F(1,2);new_s=F(2,3)
    old_a,old_b=A.substitute(s=old_s),B.substitute(s=old_s)
    new_a,new_b=A.substitute(s=new_s),B.substitute(s=new_s)
    different((old_a+old_b)**2*(4*new_a*new_b)
              -(new_a+new_b)**2*(4*old_a*old_b),'seed change alters cross-ratio')
    simple_wave=x*t
    fixed_read=simple_wave.derivative('t')+simple_wave.derivative('x')
    different(fixed_read.substitute(t=1)-fixed_read.substitute(t=0),
              'fixed-position wave quantity changes')
    require((Q*32767)//S==(Q*32768)//S,'quotient collision witness')
    require((Q*32767)%S!=(Q*32768)%S,'residual distinction')
    # Finite samples alone do not close arbitrary wave dynamics.
    vanishing=P.coerce(1)
    nodes=tuple(zi.substitute(s=F(1,2))*F(1,4) for zi in z)
    for node in nodes: vanishing*=x-node
    for node in nodes: identity(vanishing.substitute(x=node),'zero initial sample')
    evolved=evolve_polynomial(vanishing,P.coerce(0),x,t)
    different(evolved.substitute(x=nodes[-1],t=F(1,10)),
              'degree-six wave produces a different later reading')
    print('PASS 8: five failure cases detected without weakening the assumptions')
    print('RESULT: 8/8 groups passed. Algebraic identities plus declared finite execution only.')
    print('No proof-assistant certification or physical unification is claimed.')


if __name__=='__main__':
    main()