"""Exact formal checks for the Graviton geometric-matter research update.

No third-party packages. No writes to the manuscript or website.
Mixed derivatives commute; fields are smooth in a local Minkowski patch.
"""
from collections import Counter
from fractions import Fraction as Q
from itertools import product
import json


def clean(p):
    return {k:v for k,v in p.items() if v}


def derivative_exp_polynomial(p):
    # p represents exp(s*r)*sum c[s_power,r_power]*s**s_power*r**r_power.
    q=Counter()
    for (sp,rp),c in p.items():
        q[sp+1,rp]+=c
        q[sp,rp-1]+=rp*c
    return clean(q)


base={(0,-1):Q(1)}
first=derivative_exp_polynomial(base)
second=derivative_exp_polynomial(first)
radial=Counter(second)
for (sp,rp),c in first.items():
    radial[sp,rp-1]+=2*c
assert clean(radial)=={(2,-1):Q(1)} # Laplacian(exp(s*r)/r)=s*s*exp(s*r)/r, r>0.


def differentiated_h(a,b,derivs):
    # Formal derivative key: (component of A, sorted derivative indices).
    q=Counter()
    q[b,tuple(sorted((*derivs,a)))]+=1
    q[a,tuple(sorted((*derivs,b)))]+=1
    return q


def riemann_twice(a,b,c,d):
    q=Counter()
    for sign,h_a,h_b,derivs in (
        (1,a,d,(c,b)), (1,b,c,(d,a)),
        (-1,b,d,(c,a)), (-1,a,c,(d,b))
    ):
        for key,value in differentiated_h(h_a,h_b,derivs).items():
            q[key]+=sign*value
    return clean(q)


curvature=[riemann_twice(*indices) for indices in product(range(4),repeat=4)]
assert all(not component for component in curvature)

substitutions=[]
for u,L in ((Q(2),Q(1)),(Q(1,4),Q(1)),(Q(0),Q(1))):
    a=Q(3,2)*u-L
    original_residual=4*a+4*L
    corrected=4*L/a
    corrected_residual=-a*corrected+4*L
    assert original_residual==6*u
    assert corrected_residual==0
    substitutions.append({'mu_squared':str(u),'Lambda':str(L),'a':str(a),
        'original_residual':str(original_residual),'corrected_particular':str(corrected),
        'corrected_residual':str(corrected_residual)})

L=Q(1)
singular_coefficient=-Q(2,3)*L
assert 6*singular_coefficient+4*L==0 # Laplacian(r**2)=6.
assert Q(1,2)-Q(1)==-Q(1,2) # Eq9 + E=kappa*Pi differs from printed Eq15.

results={'scope':'Exact checks of displayed identities only, not verification of the full theory',
    'radial_identity':'Laplacian(exp(s*r)/r)=s^2 exp(s*r)/r for r>0',
    'formal_pure_gauge_riemann_components_checked':len(curvature),
    'nonzero_pure_gauge_riemann_components':sum(bool(x) for x in curvature),
    'constant_particular_substitution_checks':substitutions,
    'a_zero_r_squared_particular_residual':'0',
    'source_pdf_sha256':'99519795c55879827c7991a9e174a6d4317ea8f5c2df68417a3c4f6901e87033'}
print(json.dumps(results,indent=2))