generate.py

back to table · edit · history · where entries came from · files · download

3729 bytes, as of the version from 2026-10-30 06:50 (current). Recorded here, not run.

"""Radius of Lane–Emden polytropes -- numberdb.org/T613.

Run with SageMath (install numberdb once with sage -pip install numberdb):
    sage -python generate.py              # verify every stored entry
    sage -python generate.py --fill       # write entries to the existing draft
Set NUMBERDB_API_KEY for authenticated draft reads and writes. In the campaign,
NUMBERDB_KEY_FROM_STDIN=1 reads the token from stdin; NUMBERDB_PUBLISH=1 selects
entry writing. No action in this script offers, reviews or publishes a table.
Download solver.py and all files listed in `files` beside this script.

Default selection: n=0,1/2,...,9/2, for one first-zero radius per model.
See argument-selection.md for disjoint intrinsic, application and reference
roles. The solver supports these half-integer exponents, not all rational n.
The two analytic radii use Arb. Other rows use two Taylor continuations in
sqrt(theta), at 150 and 210 decimal digits, with independently varied order,
step cap and starting radius. The measured worst relative discrepancy was
4.31e-154 before decimal serialization. These are heuristic agreement checks,
not rigorous bounds for ODE truncation or error propagation.
The independent radial-coordinate implementation checks every selected model
to at least 35 digits, with exact and deliberately displaced controls.
"""
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=('0','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=32)
def compute(n,working):
    order,divisor,start=CONFIG[working]
    return str(radius(n,working,order,divisor,start))

class LaneEmdenRadius(numberdb.Generator):
    table='T613'
    parameters=('n',)
    type='R'
    digits=100
    rigour='heuristic (agreement-checked)'
    files=('generate.py','solver.py','argument-selection.md','compute.py',
           'numerical-evidence.json','independent_checks.py','independent-evidence.json',
           'readback_checks.py','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 finite half-integer index.')
        if digits>100:raise ValueError('The current precision policy was validated for 100 digits.')
        B=RealBallField(numberdb.bits(digits,losing=80))
        if n=='0':return {'number':B(6).sqrt(),'comment':r'The incompressible profile is $\theta_0(\xi)=1-\xi^2/6$, so $\xi_1(0)=\sqrt6$.'}
        if n=='1':return {'number':B.pi(),'comment':r'The linear profile is $\theta_1(\xi)=\sin(\xi)/\xi$, so $\xi_1(1)=$ HREF{Pi}[$\pi$].'}
        return numberdb.agreeing(lambda working:compute(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
    generator=LaneEmdenRadius()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1':
        # Compute before opening the streaming writer, so all ten arrive together.
        for params in generator.enumerate():
            generator.value(params, generator.digits)
        print(generator.publish(message='Add ten checked Lane-Emden radii and reproducible sources.'))
    else:
        report=generator.verify(sample=None);print(report)
        sys.exit(0 if report.ok else 1)