"""Mass-radius constants of Lane–Emden polytropes -- numberdb.org/T618.

Run with SageMath (install numberdb once with sage -pip install numberdb):
    sage -python generate.py           # recompute and verify every stored row
    sage -python generate.py --fill    # fill this existing draft
Download solver.py and the other attached files listed in `files` beside this
script. Draft reads and writes need NUMBERDB_API_KEY. The campaign wrapper
accepts NUMBERDB_KEY_FROM_STDIN=1 and NUMBERDB_PUBLISH=1 for filling entries.
Filling does not offer, review or publish the draft.

N_n = (4*pi)^(1/n)/(n+1) * mu_n^((1-n)/n) * xi_1^((n-3)/n),
so K = N_n G M^((n-1)/n) R^((3-n)/n) for P=K*rho^(1+1/n).
The index-one value is 2/pi, evaluated in Arb ball arithmetic.
Every other value integrates the normalized positive Lane–Emden solution
freshly. solver.py is the unmodified attachment from T617, itself shared by
this family; its download URL and SHA-256 are in source-manifest.json.
Runs at 150/210 decimal digits use orders 160/210, step caps 1/12 and 1/18,
and central matching radii 1/2 and 2/3. Their agreement supports only a
heuristic claim: ODE truncation, propagation and endpoint errors are not
rigorously bounded. 100 significant digits are retained.
The solver supports half-integer indices, not all rational n in the domain.
No stored numerical input is used by value(). See argument-selection.md for
the mixed argument families and the deliberately noncanonical finite bounds.
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.real_arb import RealBallField
from sage.rings.rational_field import QQ
from solver import radius

INDICES=('1/2','1','3/2','2','5/2','3','7/2','4','9/2')
CONFIG={150:(160,12,'1/2'),210:(210,18,'2/3')}

@lru_cache(maxsize=24)
def background(n,working):
    order,divisor,start=CONFIG[working]
    return radius(n,working,order,divisor,start,trace=True)

def coefficient(n,working):
    x,q,meta=background(n,working)
    F=x.parent(); nn=F(QQ(n)); one=F(1)
    return (4*F.pi())**(one/nn)/(nn+1)*q**((1-nn)/nn)*x**((nn-3)/nn)

class LaneEmdenMassRadius(numberdb.Generator):
    table='T618'
    parameters=('n',)
    type='R'
    digits=100
    rigour='heuristic (agreement-checked)'
    files=('generate.py','solver.py','compute.py','independent_checks.py',
           'radial_check.py','readback_checks.py','dry_run.py','check.py','argument-selection.md',
           'numerical-evidence.json','independent-evidence.json',
           'skill-manifest.json','skill-evidence.json','source-manifest.json')

    def enumerate(self):
        for n in INDICES:yield {'n':n}

    def value(self,params,digits):
        n=str(QQ(params['n']))
        if n not in INDICES:
            raise ValueError('Use a selected positive half-integer below five.')
        if digits>100:
            raise ValueError('This precision policy was validated for at most 100 digits.')
        if n=='1':
            B=RealBallField(numberdb.bits(digits,losing=80))
            return {'number':B(2)/B.pi(),'comment':r'The linear profile $\theta_1(\xi)=\sin(\xi)/\xi$ gives $N_1=2/\pi$.'}
        return numberdb.agreeing(lambda working:str(coefficient(n,working)),at=(150,210))

if __name__=='__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        token=sys.stdin.read().strip()
        if '=' in token and token.split('=',1)[0].replace('_','').isalnum():
            token=token.split('=',1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY']=token
    g=LaneEmdenMassRadius()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1':
        print(g.publish(message='Add nine agreement-checked polytropic mass-radius coefficients.'))
    else:
        report=g.verify(sample=None);print(report);sys.exit(0 if report.ok else 1)