"""Axial inertia of Lane–Emden polytropes -- numberdb.org/T616.

Adapted from the live T613 solver.py attachment. Adds the density moment
J=integral xi^4 theta^n dxi and H=integral xi^2 theta dxi. The identity
J=mu*xi_1^2-6H provides a separate integral check at the surface.

Only numberdb.sage initialises Sage. No symbolic power-series machinery is used.
All divisions involving computed quantities are in a named RealField.
"""
import numberdb.sage as numberdb
from sage.rings.real_mpfr import RealField
from sage.rings.rational_field import QQ
from sage.rings.integer_ring import ZZ
from math import comb

def horner(a,h):
    y=a[-1]
    for v in a[-2::-1]: y=y*h+v
    return y

def center(n,F,order,x):
    """theta=sum a[k]*xi**(2k); b=theta**n, with a[0]=b[0]=1."""
    n=F(n);a=[F(1)];b=[F(1)]
    for k in range(order):
        a.append(-b[k]/((2*k+2)*(2*k+3)))
        j=k+1
        b.append(sum(((n+1)*i-j)*a[i]*b[j-i] for i in range(1,j+1))/j)
    theta=horner(a,x*x)
    q=-x*horner([2*k*a[k] for k in range(len(a))],x*x)
    J=x**5*horner([b[k]/(2*k+5) for k in range(len(b))],x*x)
    H=x**3*horner([a[k]/(2*k+3) for k in range(len(a))],x*x)
    return theta,q,J,H

def surface_series(t,x0,q0,j0,h0,m,order):
    """x'=-2t x²/q, q'=-2t**m x⁴/q; derivatives with respect to t."""
    F=x0.parent();zero=F(0)
    x=[x0];q=[q0];js=[j0];hs=[h0];xx=[];inv=[];c=[];d=[];e=[]
    power=[F(comb(m,j))*t**(m-j) for j in range(m+1)]
    cube=[t**3,3*t*t,3*t,F(1)]
    for k in range(order):
        xx.append(sum((x[j]*x[k-j] for j in range(k+1)),zero))
        inv.append(F(1)/q0 if k==0 else -sum(q[j]*inv[k-j] for j in range(1,k+1))/q0)
        c.append(sum(xx[j]*inv[k-j] for j in range(k+1)))
        d.append(sum(xx[j]*c[k-j] for j in range(k+1)))
        e.append(sum(xx[j]*d[k-j] for j in range(k+1)))
        js.append(-2*sum(power[j]*e[k-j] for j in range(min(m,k)+1))/(k+1))
        hs.append(-2*sum(cube[j]*d[k-j] for j in range(min(3,k)+1))/(k+1))
        x.append(-2*(t*c[k]+(c[k-1] if k else zero))/(k+1))
        q.append(-2*sum(power[j]*d[k-j] for j in range(min(m,k)+1))/(k+1))
    return x,q,js,hs

def inertia(n,working=150,order=160,divisor=12,start='1/2',trace=False):
    """Return a fixed-precision approximation; there is no certified ODE bound.

    Tail-based steps use the final eight Taylor coefficients. Their predicted
    contribution is <10**(-(working-20)); the cap also changes with divisor.
    Precision, order, initial matching radius and step cap are varied in checks.
    """
    n=QQ(n);m=int(2*n+1)
    if n<0 or n>=5 or 2*n not in ZZ:
        raise ValueError('This solver supports integer/half-integer 0 <= n < 5')
    F=RealField(numberdb.bits(working));x=F(QQ(start))
    theta,q,J,H=center(n,F,2*working,x)
    t=theta.sqrt();count=0;tol=F(10)**(-(working-20));max_tail=F(0)
    while t>0:
        a,b,c,d=surface_series(t,x,q,J,H,m,order)
        # Center singularity at t=1; possible distant complex singularities
        # are additionally detected by the tail estimate, not assumed absent.
        cap=min(F(1)/divisor,(1-t)/3,t)
        for arr,scale in [(a,max(F(1),abs(x))),(b,max(F(1),abs(q))),(c,max(F(1),abs(J))),(d,max(F(1),abs(H)))]:
            for k in range(order-7,order+1):
                if arr[k]:
                    cap=min(cap,(tol*scale/abs(arr[k]))**(F(1)/k)/2)
        h=-cap
        xn=horner(a,h);qn=horner(b,h)
        tail=max(abs(a[-1]*h**order)/max(F(1),abs(xn)),abs(b[-1]*h**order)/max(F(1),abs(qn)))
        max_tail=max(max_tail,tail)
        assert xn>x and qn>0 and qn>=q*(1-F(10)**(-working+10))
        J,H=horner(c,h),horner(d,h)
        x,q=xn,qn;t+=h;count+=1
        if count>3000:raise RuntimeError('too many continuation steps')
    meta={'steps':count,'last_term_max':str(max_tail),'working':working,'order':order,'divisor':divisor,'start':start}
    value=2*J/(3*q*x*x)
    alternate=F(2)/3-4*H/(q*x*x)
    discrepancy=abs(value-alternate)/abs(value)
    assert discrepancy<F(10)**(-working+30), ('integral identity',n,discrepancy)
    meta.update(radius=str(x),mass=str(q),density_moment=str(J),enthalpy_moment=str(H),identity_relative_error=str(discrepancy))
    return (value,meta) if trace else value