"""
A small exact polynomial normalizer, not a proof-assistant kernel.
Coefficients are Rational{BigInt}. Two explicitly justified relations are used:
    phi^2 = phi + 1,   c^2 = 1 - s^2.
Each monomial reduction substitutes one of these identities, so a zero
normal form is a sufficient certificate of a polynomial identity for all
real assignments satisfying the relations. No numerical substitution is used.
Nonzero denominator and positivity conditions are supplied separately in P01-P12.
"""
const POLY_VARS = (:a,:phi,:pi,:s,:c,:z,:w,:x,:k,:l,:alpha,:b,:u,:v,
                   :d,:f,:r,:h,:e,:A,:B,:t,:V,:An,:Bn,:ell)
const Monomial = NTuple{length(POLY_VARS),Int}
const ZERO_MONOMIAL = ntuple(_->0,length(POLY_VARS))
const PHI_INDEX = 2
const SIN_INDEX = 4
const COS_INDEX = 5

struct Poly
    terms::Dict{Monomial,Rat}
end
Poly(x::Union{Integer,Rational}) = x == 0 ? Poly(Dict{Monomial,Rat}()) :
    Poly(Dict(ZERO_MONOMIAL=>rational(x)))
function variable(name::Symbol)
    i=findfirst(==(name),POLY_VARS)
    isnothing(i) && throw(ArgumentError("unknown polynomial variable: $name"))
    Poly(Dict(ntuple(j->Int(i==j),length(POLY_VARS))=>rational(1)))
end
function addterm!(out::Dict{Monomial,Rat},m::Monomial,c::Rat)
    c == 0 && return out
    v=get(out,m,rational(0))+c
    v == 0 ? delete!(out,m) : (out[m]=v)
    out
end
function reduced_term!(out::Dict{Monomial,Rat},m::Monomial,c::Rat)
    c == 0 && return out
    if m[PHI_INDEX] >= 2
        n=Base.setindex(m,m[PHI_INDEX]-2,PHI_INDEX)
        reduced_term!(out,n,c)
        reduced_term!(out,Base.setindex(n,n[PHI_INDEX]+1,PHI_INDEX),c)
    elseif m[COS_INDEX] >= 2
        n=Base.setindex(m,m[COS_INDEX]-2,COS_INDEX)
        reduced_term!(out,n,c)
        reduced_term!(out,Base.setindex(n,n[SIN_INDEX]+2,SIN_INDEX),-c)
    else
        addterm!(out,m,c)
    end
    out
end
Base.iszero(p::Poly) = isempty(p.terms)
Base.zero(::Type{Poly}) = Poly(0)
Base.zero(::Poly) = Poly(0)
Base.one(::Type{Poly}) = Poly(1)
Base.one(::Poly) = Poly(1)
Base.:-(p::Poly) = Poly(Dict(m=>-c for (m,c) in p.terms))
function Base.:+(p::Poly,q::Poly)
    out=copy(p.terms)
    for (m,c) in q.terms; addterm!(out,m,c); end
    Poly(out)
end
Base.:-(p::Poly,q::Poly) = p+(-q)
function Base.:*(p::Poly,q::Poly)
    out=Dict{Monomial,Rat}()
    for (m,c) in p.terms, (n,d) in q.terms
        reduced_term!(out,ntuple(i->m[i]+n[i],length(POLY_VARS)),c*d)
    end
    Poly(out)
end
for op in (:+,:-,:*)
    @eval Base.$op(p::Poly,x::Union{Integer,Rational}) = Base.$op(p,Poly(x))
    @eval Base.$op(x::Union{Integer,Rational},p::Poly) = Base.$op(Poly(x),p)
end
Base.:/(p::Poly,x::Union{Integer,Rational}) = p*inv(rational(x))
function Base.:^(p::Poly,n::Integer)
    n >= 0 || throw(ArgumentError("polynomial exponent must be nonnegative"))
    result=Poly(1); base=p; k=n
    while k > 0
        isodd(k) && (result=result*base)
        k=div(k,2)
        k > 0 && (base=base*base)
    end
    result
end
Base.:(==)(p::Poly,q::Poly) = iszero(p-q)

"""Canonical readable text; exponent order is fixed and deterministic."""
function Base.show(io::IO,p::Poly)
    if iszero(p); print(io,"0"); return; end
    texts=String[]
    for m in sort(collect(keys(p.terms)))
        factors=String[string(p.terms[m])]
        for (i,e) in enumerate(m)
            e == 0 && continue
            push!(factors,e==1 ? string(POLY_VARS[i]) : "$(POLY_VARS[i])^$e")
        end
        push!(texts,join(factors,"*"))
    end
    print(io,join(texts," + "))
end

function polynomial_certificates()
    a,phi,pi,s,c,z,w,x,k,l,alpha,b,u,v,d,f,r,h,e,A,B,t,V,An,Bn,ell = variable.(POLY_VARS)
    out=Pair{String,Poly}[]
    put(name,p)=push!(out,name=>p)
    eps=pi/648000; delta=eps/2
    put("phi_sq",phi^2-phi-1)
    put("sc_identity",s^2+c^2-1)
    put("one_arcsecond_cleared",eps*180*3600-pi)
    put("shares_sum_cleared",phi+pi-(phi+pi))
    put("mu_eq_one_sub_cleared",pi-((phi+pi)-phi))
    put("weighted_mean_cleared",(a*phi*x+a*pi)-(a*(phi+pi)+a*phi*(x-1)))
    Tk(k,x)=1+k*(x-1)
    H(k,x)=alpha+k*(x-alpha)
    put("displacement_conjugacy",Tk(k,x)-1-k*(x-1))
    put("affine_iterate_base",Tk(1,x)-x)
    put("affine_iterate_step",Tk(k,Tk(l,x))-Tk(k*l,x))
    put("homothety_comp",H(k,H(l,z))-H(k*l,z))
    put("homothety_one",H(1,z)-z)
    put("observation_commutes_after_linearity",1+k*ell-Tk(k,1+ell))
    put("arbitrary_intercept_defect",b+k*ell-Tk(k,b+ell)-(1-k)*(b-1))
    put("grid_interface_cleared",k*(h*d)-(k*d)*h)
    put("mark_distance_squared",(a*c-a*c)^2+(a*z-a*w)^2-a^2*(z-w)^2)
    put("inner_extension_offset",(1+s)-s-1)
    put("outer_extension_offset",(phi^2+s)-(1+s)-phi)
    put("negative_outer_extension_offset",-(phi^2+s)+(1+s)+phi)
    put("total_extension_offset",phi^2+s-s-phi^2)
    put("chord_offset",s-(-s)-2*s)
    put("sequential_upper_inner",a*s+a-a*(s+1))
    put("sequential_upper_outer",a*(s+1)+a*phi-a*(s+phi^2))
    put("sequential_lower_inner",-a*s-a+a*(s+1))
    put("sequential_lower_outer",-a*(s+1)-a*phi+a*(s+phi^2))
    put("G_eq",(a*phi+a*phi)/2-a*phi)
    put("C_eq",a*(3*pi/2-pi/2)-a*pi)
    put("W_eq",a*phi+a*pi-a*(phi+pi))
    put("geometric_share_G_cleared",a*phi*(phi+pi)-phi*(a*phi+a*pi))
    put("geometric_share_C_cleared",a*pi*(phi+pi)-pi*(a*phi+a*pi))
    put("outline_length",a*(2*pi-delta-delta)+2*a*phi^2-(a*(2*pi-eps)+2*a*phi^2))
    put("reference_not_outline_positive_factor",a*(2*pi-eps)+2*a*phi^2-a*(phi+pi)-a*(pi-eps+phi+2))
    heights=[-(phi^2+s),-(1+s),-s,s,1+s,phi^2+s]
    put("marks_collinear_x",a*c-a*c)
    put("seeds_distinct_difference",a*s-(-a*s)-2*a*s)
    # Seed coefficients are (x/(ac) +/- y/(as))/2. Clear 2acs.
    put("seed_linear_combination_re_cleared",(x*s+v*c)*(a*c)+(x*s-v*c)*(a*c)-2*a*c*s*x)
    put("seed_linear_combination_im_cleared",(x*s+v*c)*(a*s)-(x*s-v*c)*(a*s)-2*a*c*s*v)
    put("seed_minor",(a*c)*(-a*s)-(a*c)*(a*s)+2*a^2*c*s)
    put("moment_11",sum([(a*c)^2 for _ in heights])-6*(a*c)^2)
    put("moment_12",sum([a*c*a*zi for zi in heights]))
    put("moment_22",sum([(a*zi)^2 for zi in heights])-2*a^2*(s^2+(1+s)^2+(phi^2+s)^2))
    for i in 1:6, j in 1:6
        put("gram_entry_$(i)_$(j)",(a*c)^2+(a*heights[i])*(a*heights[j])-a^2*(c^2+heights[i]*heights[j]))
    end
    put("mark_dot",(a*c)^2+(a*z)*(a*w)-a^2*(c^2+z*w))
    put("mark_radius_squared",(a*c)^2+(a*z)^2-a^2*(c^2+z^2))
    put("seed_angle_cosine",(c^2-s^2)-(2*c^2-1))
    put("triangle_area_signed_determinant",a*c*a*w-a*z*a*c-a^2*c*(w-z))
    put("circle_speed_squared",(-a*s)^2+(a*c)^2-a^2)
    put("arc_antiderivative_endpoints",a*v-a*u-a*(v-u))
    put("mark_scale_re",k*(a*c)-(k*a)*c)
    put("mark_scale_im",k*(a*z)-(k*a)*z)
    put("circle_scale_re",k*(a*c)-(k*a)*c)
    put("circle_scale_im",k*(a*s)-(k*a)*s)
    put("placed_homothety_after_linearity",alpha+k*((alpha+ell)-alpha)-(alpha+k*ell))
    put("readout_scale_G",(k*a)*phi-k*(a*phi))
    put("readout_scale_C",(k*a)*pi-k*(a*pi))
    put("farey_determinant_identity",20706*(12343*k-4196*d)+12343*(7039*d-20706*k)-d)
    put("rational_residual_cleared",A*d-B*f-(A*d-B*f))
    put("real_residual_after_euclidean",k*d-f-((q*d-f)+(k-q)*d))
    put("scalar_roundtrip",h*e+(z-h*e)-z)
    put("scalar_canonical_residual",h*e+r-h*e-r)
    put("ledger_one_step_cleared",(B*t+r)-(B*t+r))
    # An = A^n, Bn = B^n, V is the tail polynomial. These are algebraically
    # independent here; the induction substitutes their meanings afterwards.
    put("ledger_weighted_induction",B*(Bn*t+V)+An*r-(Bn*B*t+(An*r+B*V)))
    put("ledger_step_distribute",An*(B*t+r)-(B*(An*t)+An*r))
    put("point_reconstruction_substitution",alpha+ell*(h*d)-(alpha+ell*(h*d)))
    out
end