# OMEGA Geometry -> Center-relative Ellipse -> Hilbert-space Operator
# Version 1.0, 2026-10-07. Julia 1.9+ target; standard libraries only.
# EXECUTION STATUS AT DELIVERY: Julia runtime not available; NOT RUN in Julia.
# See START_HERE.html for proofs, source attribution, preservation contracts,
# the independent Python execution record, and the translated-tail correction.
# Mathematical equalities in the written proof are exact. Float64 checks below
# are numerical conformance checks, not interval certificates or a Lean proof.
# The four Q(sqrt(5))[s] checks use exact rational arithmetic.
# Source: https://omega-naught.grok.me/papers/OMEGA_Paper_VII_v1.6.0.pdf
# Source DOI: 10.5281/zenodo.22233573, definitions 3.1-3.3, equations 39-46.
# The moment ellipse and direct-sum extensions are NEW DERIVATIONS here;
# they are not claimed to have been stated in Paper VII.

module OmegaGeometryBridge
using LinearAlgebra

export kernel_record, arc_point, foot_point, place, moment_ellipse, ellipse_affine_chart,
       ellipse_operator, numerical_boundary, ellipse_witness,
       block_truncation, truncation_error, recover_family,
       matrix_polynomial, scalar_polynomial, run_checks

require(ok::Bool, message::AbstractString) = ok || throw(ArgumentError(message))
function near(x, y; atol=2e-12, rtol=2e-11)
    return norm(x-y) <= atol + rtol * max(norm(x), norm(y))
end
finitecomplex(z::Number) = isfinite(real(z)) && isfinite(imag(z))

"""Canonical Paper VII geometry by default. Other epsilon values define an
explicit auxiliary family, not a change to the published one-arcsecond seed.
X is in the intrinsic basis (v,u), not canonical world coordinates (x,y,z).
"""
function kernel_record(a::Real=1.0; epsilon::Real=pi/648000)
    aa, ee = Float64(a), Float64(epsilon)
    require(isfinite(aa) && aa > 0, "a must be finite and positive")
    require(isfinite(ee) && 0 < ee < pi, "require 0 < epsilon < pi")
    phi = (1 + sqrt(5.0))/2
    delta = ee/2
    s, c = sin(delta), cos(delta)
    require(s > 0 && c > 0, "seed underflow or degenerate geometry")
    z = [-phi^2-s, -1-s, -s, s, 1+s, phi^2+s]
    X = aa .* vcat(reshape(z, 1, 6), fill(c, 1, 6))
    require(all(isfinite, X), "coordinate overflow")
    G = transpose(X)*X
    require(all(isfinite, G), "Gram overflow")
    bearings = atan.(z, c)
    angles = [abs(bearings[j]-bearings[i]) for i in 1:6, j in 1:6]
    world = vcat(reshape(aa .* z, 1, 6), fill(-aa*c,1,6), zeros(1,6))
    return (a=aa, epsilon=ee, delta=delta, phi=phi, s=s, c=c,
            z=z, X=X, G=G, bearings=bearings, angles=angles, world=world)
end

function arc_point(k, t::Real)
    tt = Float64(t)
    require(isfinite(tt) && k.delta <= tt <= 2pi-k.delta,
            "t must belong to the retained major arc")
    return k.a .* [sin(tt), cos(tt)]
end
function foot_point(k, sign::Integer, ell::Real)
    ll = Float64(ell)
    require(sign == -1 || sign == 1, "sign must be -1 or 1")
    require(isfinite(ll) && 0 <= ll <= k.phi^2,
            "ell must be in [0,phi^2]")
    return k.a .* [sign*(k.s+ll), k.c]
end

"""Place intrinsic columns with a numerically orthonormal m-by-2 frame.
The analytic theorem applies to exact isometries, in arbitrary Hilbert space.
"""
function place(Y::AbstractMatrix, E::AbstractMatrix, alpha::AbstractVector)
    require(size(Y,1)==2 && size(E,2)==2 && size(E,1)==length(alpha),
            "incompatible coordinate/frame dimensions")
    require(all(isfinite,Y) && all(isfinite,E) && all(isfinite,alpha),
            "nonfinite placement input")
    require(near(transpose(E)*E, Matrix{Float64}(I,2,2)),
            "frame columns must be orthonormal within numerical tolerance")
    return reshape(alpha, :, 1) .+ E*Y
end

"""E_Omega = (X/sqrt(6))(closed unit ball in R^6).
Equal weights and the circle center as origin are explicit choices.
This ellipse is NOT the arc-and-feet curve and NOT centroid covariance.
"""
function moment_ellipse(k)
    R = k.a * norm(k.z)/sqrt(6.0)
    r = k.a * k.c
    require(isfinite(R) && isfinite(r) && R > r > 0,
            "unrepresentable/non-elliptic semiaxes")
    return (R=R, r=r, L=k.X/sqrt(6.0), M=(k.X*transpose(k.X))/6)
end

"""Map the complete intrinsic kernel through D. Its original circle maps
onto the moment ellipse boundary. Ordinary Euclidean angles generally change;
H is the target metric that recovers the source Euclidean inner products.
"""
function ellipse_affine_chart(k)
    e=moment_ellipse(k)
    D=Diagonal([e.R/k.a,e.r/k.a])
    H=Diagonal([(k.a/e.R)^2,(k.a/e.r)^2])
    return (D=D,H=H,marks=D*k.X)
end

"""Construct zeta*I + exp(i*angle)*[d 2r; 0 -d].
Its numerical range is the FILLED translated/rotated ellipse.
"""
function ellipse_operator(R::Real, r::Real; center::Number=0+0im,
                          angle::Real=0.0)
    RR, rr, th, zz = Float64(R), Float64(r), Float64(angle), ComplexF64(center)
    require(isfinite(RR) && isfinite(rr) && RR >= rr > 0,
            "require finite R >= r > 0")
    require(isfinite(th) && finitecomplex(zz), "nonfinite center/angle")
    # Factored expression avoids cancellation in R^2-r^2 near a circle.
    dd = sqrt(RR-rr)*sqrt(RR+rr)
    require(isfinite(dd) && isfinite(2rr), "operator overflow")
    T0 = ComplexF64[dd 2rr; 0 -dd]
    T = zz .* Matrix{ComplexF64}(I,2,2) .+ cis(th).*T0
    require(all(finitecomplex,T), "translated operator overflow")
    return T
end

function numerical_boundary(T::AbstractMatrix, theta::Real)
    require(size(T,1)==size(T,2) && size(T,1)>0, "T must be square")
    require(all(finitecomplex,T) && isfinite(theta), "nonfinite boundary input")
    q = cis(-Float64(theta))
    H = Hermitian((q*T + conj(q)*adjoint(T))/2)
    F = eigen(H)
    v = F.vectors[:,end]
    return (value=dot(v,T*v), support=F.values[end], vector=v)
end

"""Numerical unit-vector witness for a point of a centered ellipse.
A roundoff band of 256*eps(Float64) is allowed in normalized membership.
An exact membership decision is not claimed.
"""
function ellipse_witness(R::Real, r::Real, w::Number)
    RR, rr, ww = Float64(R), Float64(r), ComplexF64(w)
    require(isfinite(RR) && isfinite(rr) && RR >= rr > 0 && finitecomplex(ww),
            "invalid ellipse/witness input")
    dd = sqrt(RR-rr)*sqrt(RR+rr)
    xx, yy = real(ww), imag(ww)
    t = (xx/RR)^2 + (yy/rr)^2
    require(t <= 1 + 256eps(Float64), "point outside ellipse")
    q0 = [rr*xx/RR^2, yy/rr, dd*xx/RR^2]
    n = [dd/RR, 0.0, -rr/RR]
    X,Y,Z = q0 + sqrt(max(0.0,1-t))*n
    if Z >= 0
        v1 = sqrt((1+Z)/2)
        v = ComplexF64[v1, (X+im*Y)/(2v1)]
    else
        v2 = sqrt((1-Z)/2)
        v = ComplexF64[(X-im*Y)/(2v2), v2]
    end
    return v
end

"""First N blocks of zeta*I + exp(i*angle)*direct_sum(2^(1-n)*T0).
The omitted infinite tail of the TRUNCATED operator is zeta*I, NOT zero
unless zeta=0. This finite matrix does not numerically simulate infinity.
"""
function block_truncation(R::Real, r::Real, N::Integer;
                          center::Number=0+0im, angle::Real=0.0)
    require(1 <= N <= 512, "require 1 <= N <= 512 for dense truncation")
    T0 = ellipse_operator(R,r)
    zz, th = ComplexF64(center), Float64(angle)
    require(finitecomplex(zz) && isfinite(th), "nonfinite center/angle")
    B = zz .* Matrix{ComplexF64}(I,2N,2N)
    for n in 1:N
        ii = (2n-1):(2n)
        B[ii,ii] .+= (cis(th)*ldexp(1.0,1-n)) .* T0
    end
    return B
end
function truncation_error(R::Real,r::Real,N::Integer)
    require(isfinite(R) && isfinite(r) && R >= r > 0 && N >= 0,
            "invalid tail-bound parameters")
    # Float64 evaluation of the exact formula (R+r)*2^(-N).
    return ldexp(Float64(R)+Float64(r), -N)
end

"""Recover s and a within the declared symmetric golden-step family only.
Requires squared aspect ratio > phi^2. No arbitrary-shape inverse is claimed.
"""
function recover_family(R::Real, r::Real)
    RR, rr = Float64(R), Float64(r)
    phi = (1+sqrt(5.0))/2
    require(isfinite(RR) && isfinite(rr) && RR >= rr > 0,
            "invalid semiaxes")
    q = (RR/rr)^2
    require(isfinite(q) && q > phi^2, "outside declared family's aspect ratios")
    h = 2*(phi+2)/3
    s = 2*(q-phi^2)/(h+sqrt(h*h+4*(q+1)*(q-phi^2)))
    require(0 < s < 1, "unrepresentable recovered parameter")
    a = rr/sqrt((1-s)*(1+s))
    return (s=s, epsilon=2asin(s), a=a)
end

function scalar_polynomial(coeffs::AbstractVector, z::Number)
    require(!isempty(coeffs), "empty coefficient list")
    m = size(first(coeffs),1)
    require(m > 0 && all(C -> size(C)==(m,m) && all(finitecomplex,C),coeffs),
            "coefficient matrices must be finite, square and equally sized")
    require(finitecomplex(z), "nonfinite polynomial argument")
    V = zeros(ComplexF64,m,m)
    for C in reverse(coeffs)
        V = V*z + C
    end
    return V
end
function matrix_polynomial(coeffs::AbstractVector, T::AbstractMatrix)
    require(size(T,1)==size(T,2) && size(T,1)>0 && all(finitecomplex,T),
            "invalid operator argument")
    scalar_polynomial(coeffs,0) # validate coefficient shapes
    m,n = size(first(coeffs),1),size(T,1)
    V = zeros(ComplexF64,m*n,m*n)
    power = Matrix{ComplexF64}(I,n,n)
    for C in coeffs
        V += kron(C,power)
        power = power*T
    end
    return V
end

# Exact certificates in Q(sqrt(5))[s], independent of trigonometric evaluation.
struct Q5
    a::Rational{BigInt}
    b::Rational{BigInt}
end
q5(a::Integer=0,b::Integer=0,d::Integer=1) = Q5(BigInt(a)//BigInt(d),BigInt(b)//BigInt(d))
Base.:+(x::Q5,y::Q5) = Q5(x.a+y.a,x.b+y.b)
Base.:-(x::Q5,y::Q5) = Q5(x.a-y.a,x.b-y.b)
Base.:*(x::Q5,y::Q5) = Q5(x.a*y.a+5*x.b*y.b,x.a*y.b+x.b*y.a)
Base.:(==)(x::Q5,y::Q5) = x.a==y.a && x.b==y.b
function padd(p::Vector{Q5},q::Vector{Q5})
    r=fill(q5(),max(length(p),length(q)))
    for i in eachindex(p); r[i]=r[i]+p[i]; end
    for i in eachindex(q); r[i]=r[i]+q[i]; end
    return r
end
function pmul(p::Vector{Q5},q::Vector{Q5})
    r=fill(q5(),length(p)+length(q)-1)
    for i in eachindex(p),j in eachindex(q); r[i+j-1]=r[i+j-1]+p[i]*q[j]; end
    return r
end

function run_checks()
    names = String[]
    function group(f::Function,name::String)
        f()
        push!(names,name)
        println("PASS ",lpad(length(names),2,'0')," ",name)
    end
    function rejects(f::Function)
        try
            f()
        catch err
            err isa ArgumentError && return true
            rethrow()
        end
        return false
    end
    println("OMEGA geometry/operator bridge v1.0")
    println("Julia ",VERSION)
    println("Float64 tests are conformance checks, not formal proofs.")
    p=q5(1,1,2); one=q5(1); zero=q5()
    group("exact golden identity") do
        require(p*p==p+one,"phi^2")
    end
    group("exact fourth-power identity") do
        require(one+(p*p)*(p*p)==q5(3)*(p*p),"phi^4")
    end
    squares=padd(padd(pmul([zero,one],[zero,one]),
                          pmul([one,one],[one,one])),
                          pmul([p*p,one],[p*p,one]))
    Rsq=[q5(1,0,3)*x for x in squares]
    group("exact moment semiaxis polynomial") do
        require(Rsq==[p*p,q5(2,0,3)*(p+q5(2)),one],"R^2 polynomial")
    end
    group("exact focal polynomial") do
        require(padd(Rsq,[q5(-1),zero,one])==
                [p,q5(2,0,3)*(p+q5(2)),q5(2)],"d^2 polynomial")
    end
    k=kernel_record(); e=moment_ellipse(k); T=ellipse_operator(e.R,e.r)
    group("published unit-node transcription") do
        xs=[-2.6180364128183005,-1.0000024240684056,
            -2.424068405545306e-6,2.424068405545306e-6,
            1.0000024240684056,2.6180364128183005]
        require(near(k.world[1,:],xs;atol=2e-15,rtol=2e-15),"source x")
        require(near(k.world[2,:],fill(-0.9999999999970619,6);
                     atol=2e-15,rtol=2e-15),"source y")
    end
    group("seed, feet, retained arc and five gaps") do
        require(near(norm(k.X[:,3]),1) && near(norm(k.X[:,4]),1),"seed radii")
        require(near(diff(k.z),[k.phi,1,2k.s,1,k.phi]),"gaps")
        require(near(arc_point(k,k.delta),k.X[:,4]),"arc b")
        require(near(arc_point(k,2pi-k.delta),k.X[:,3]),"arc d")
        require(near(foot_point(k,-1,k.phi^2),k.X[:,1]),"left foot")
        require(near(foot_point(k,1,k.phi^2),k.X[:,6]),"right foot")
    end
    group("center-relative Gram and all fifteen angle cosines") do
        G=k.c^2*ones(6,6)+k.z*transpose(k.z)
        require(near(G,k.G),"Gram factorization")
        for i in 1:6,j in (i+1):6
            require(near(cos(k.angles[i,j]),G[i,j]/sqrt(G[i,i]*G[j,j])),"angle")
        end
    end
    group("angle identities without small-angle acos cancellation") do
        require(near(k.bearings[4],k.delta;atol=2e-18,rtol=2e-12),"theta0")
        require(near(k.bearings[5],pi/4+k.epsilon/4),"theta1")
        require(near(k.angles[3,4],k.epsilon;atol=2e-18,rtol=2e-12),"seed angle")
    end
    group("rank two at circle center; rank one at mark centroid") do
        require(rank(k.X)==2,"center rank")
        Xc=k.X .- sum(k.X;dims=2)/6
        require(rank(Xc)==1,"centroid rank")
    end
    group("finite isometric placement preserves full-curve samples") do
        E=zeros(9,2); E[1,1]=1/sqrt(2);E[3,1]=1/sqrt(2)
        E[2,2]=1/sqrt(2);E[9,2]=-1/sqrt(2)
        alpha=collect(0.0:8.0)/10
        Y=hcat(k.X,arc_point(k,pi),foot_point(k,-1,0.3))
        P=place(Y,E,alpha); C=P .- reshape(alpha,:,1)
        require(near(transpose(C)*C,transpose(Y)*Y),"isometric Gram")
    end
    group("ellipse matrix, singular values and constructive ball preimages") do
        require(near(e.M,Diagonal([e.R^2,e.r^2])),"moment matrix")
        require(near(svdvals(e.L),[e.R,e.r]),"semiaxes")
        for t in range(0,2pi;length=65)
            y=[e.R*cos(t),e.r*sin(t)]
            q=transpose(e.L)*(e.M\y)
            require(near(e.L*q,y) && near(norm(q),1),"ball preimage")
        end
    end
    group("whole-kernel affine ellipse chart and angle-recovery metric") do
        chart=ellipse_affine_chart(k)
        require(near(transpose(chart.marks)*chart.H*chart.marks,k.G),"metric Gram")
        for t in range(k.delta,2pi-k.delta;length=129)
            y=chart.D*arc_point(k,t)
            require(near((y[1]/e.R)^2+(y[2]/e.r)^2,1),"elliptical arc")
        end
        x=chart.D*[1.0,1.0];y=chart.D*[1.0,-1.0]
        require(abs(dot(x,y))>0.1,"Euclidean-angle distortion counterexample")
        require(near(dot(x,chart.H*y),0),"metric angle recovery")
    end
    group("scale covariance") do
        for a in (0.1,1.0,3.25,10.0)
            ka=kernel_record(a);ea=moment_ellipse(ka)
            require(near(ea.R,a*e.R) && near(ea.r,a*e.r),"axes scaling")
            require(near(ka.G,a^2*k.G),"Gram scaling")
            require(near(ellipse_operator(ea.R,ea.r),a*T),"operator scaling")
        end
    end
    group("declared-family shape and scale recovery") do
        for epsilon in (pi/648000,0.2,1.2), a in (0.3,1.0,4.0)
            kk=kernel_record(a;epsilon=epsilon);ee=moment_ellipse(kk)
            rr=recover_family(ee.R,ee.r)
            require(near(rr.s,kk.s;atol=3e-15,rtol=2e-12),"s recovery")
            require(near(rr.a,a),"a recovery")
        end
    end
    group("all-point numerical-range witnesses on a sample grid") do
        for (R,r) in ((e.R,e.r),(3.0,2.0),(1.0,1.0))
            A=ellipse_operator(R,r)
            for radius in (0.0,0.2,0.8,1.0),t in range(0,2pi;length=129)
                w=radius*(R*cos(t)+im*r*sin(t));v=ellipse_witness(R,r,w)
                require(near(norm(v),1) && near(dot(v,A*v),w),"unit witness")
            end
        end
    end
    group("support and boundary including translated ellipses") do
        for (R,r,zeta,angle) in ((e.R,e.r,0+0im,0.0),
                                (5.0,1.0,2-3im,0.37),(1.0,1.0,-1+2im,1.1))
            A=ellipse_operator(R,r;center=zeta,angle=angle)
            for t in range(0,2pi;length=721)
                b=numerical_boundary(A,t);u=t-angle
                h=hypot(R*cos(u),r*sin(u))
                w=zeta+cis(angle)*(R^2*cos(u)/h+im*r^2*sin(u)/h)
                require(near(b.value,w),"boundary point")
                require(near(b.support,real(cis(-t)*zeta)+h),"support")
            end
        end
    end
    group("squared operator, norm and explicit similarity") do
        d=sqrt(e.R-e.r)*sqrt(e.R+e.r)
        S=ComplexF64[1 -e.r/e.R;0 d/e.R]
        require(near(T*T,d^2*Matrix{ComplexF64}(I,2,2)),"T^2")
        require(near(opnorm(T),e.R+e.r),"operator norm")
        require(near(T*S,S*Diagonal([d,-d])),"similarity")
        require(near(cond(S),sqrt((e.R+e.r)/(e.R-e.r))),"condition number")
    end
    group("translated zero-tail regression and corrected center tail") do
        A=ellipse_operator(1,1;center=3)
        bad=zeros(ComplexF64,3,3);bad[1:2,1:2]=A
        good=copy(bad);good[3,3]=3
        require(near(numerical_boundary(bad,pi).support,0),"bad contains zero")
        require(near(numerical_boundary(good,pi).support,-2),"good left edge")
    end
    group("direct-sum supports and exact-form truncation errors") do
        for N in (1,2,5,9)
            B=block_truncation(e.R,e.r,N;center=2-3im,angle=0.37)
            for t in range(0,2pi;length=65)
                expected=real(cis(-t)*(2-3im))+hypot(e.R*cos(t-0.37),e.r*sin(t-0.37))
                require(near(numerical_boundary(B,t).support,expected),"block support")
            end
            Bn=block_truncation(e.R,e.r,N+1)
            B0=zeros(ComplexF64,2(N+1),2(N+1));B0[1:2N,1:2N]=block_truncation(e.R,e.r,N)
            require(near(opnorm(Bn-B0),truncation_error(e.R,e.r,N)),"tail norm")
        end
    end
    group("matrix-polynomial similarity bound, no Crouzeix oracle") do
        coeffs=[ComplexF64[1 0;0 2],ComplexF64[0.2 im/3;-0.1 0.4],
                ComplexF64[0.1 -0.2;0.3 0.05]]
        d=sqrt(e.R-e.r)*sqrt(e.R+e.r)
        kappa=sqrt((e.R+e.r)/(e.R-e.r))
        for N in (1,3,6)
            A=block_truncation(e.R,e.r,N)
            bound=maximum(opnorm(scalar_polynomial(coeffs,sgn*ldexp(d,1-n)))
                          for n in 1:N for sgn in (-1,1))
            require(opnorm(matrix_polynomial(coeffs,A)) <= kappa*bound+2e-11,
                    "explicit similarity inequality")
        end
    end
    group("invalid-input rejection") do
        require(rejects(() -> kernel_record(-1)),"negative scale")
        require(rejects(() -> kernel_record(1;epsilon=0)),"zero seed")
        require(rejects(() -> ellipse_operator(1,2)),"reversed axes")
        require(rejects(() -> ellipse_operator(1,1;center=Inf)),"infinite center")
        require(rejects(() -> ellipse_operator(1,1;angle=NaN)),"NaN angle")
        require(rejects(() -> ellipse_witness(1,1,2)),"outside point")
        require(rejects(() -> block_truncation(1,1,0)),"zero blocks")
    end
    println("GROUPS ",length(names)," PASS")
    println("CANONICAL R=",e.R," r=",e.r," norm(T)=",opnorm(T))
    println("No execution of an infinite object or proof-assistant build is claimed.")
    return names
end
end # module

if abspath(PROGRAM_FILE) == abspath(@__FILE__)
    OmegaGeometryBridge.run_checks()
end