"""High-order Taylor integration of finite Lane–Emden models.

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)
    return theta,q

def surface_series(t,x0,q0,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];xx=[];inv=[];c=[];d=[]
    power=[F(comb(m,j))*t**(m-j) for j in range(m+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)))
        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

def radius(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=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=surface_series(t,x,q,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)))]:
            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))
        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}
    return (x,q,meta) if trace else x