"""Full spectral zeta values of homogeneous oscillators -- numberdb.org/T541.

Run with SageMath and numberdb installed:
    sage -python generate.py              # verify the private or public table
    sage -python generate.py --fill       # fill the draft, with NUMBERDB_API_KEY

Download accelerate.py alongside this file. This computes the full-line traces
of -d^2/dx^2+|x|^N with unit coefficients, zero shift, and multiplicity one.
For N>=3, sum Watkins (2012), equations (2.5a,b), at sigma=2/(N+2),
lambda=+/-1/2. The factor 1/4 in (2.5a) is essential. The prefactor in
his introductory equation (1.3) is inconsistent with (2.5a) and is not used.

For order two, the generalized hypergeometric sum is accelerated with a
Gamma-ratio asymptotic tail. We compare 150 decimal digits, cutoff 160 and
140 tail corrections against 200 digits, cutoff 320 and 180 corrections.
These choices were measured to agree to over 125 digits on this selection.
The result is heuristic agreement, not an enclosure of the asymptotic remainder.
N=1 uses Airy initial values; N=2 uses the exact odd-integer spectrum.
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
import mpmath as mp
from accelerate import hyper_tail

SETTINGS = {150: (160, 140), 200: (320, 180)}

@lru_cache(None)
def boundary(N, n, sign, working):
    with mp.workdps(working):
        s=mp.mpf(2)/(N+2)
        lam=mp.mpf(sign)/2
        half=mp.mpf(1)/2
        if n==1:
            x=s**(2-2*s)*mp.gamma(s*(1+lam))*mp.gamma(s)*mp.gamma(half-s)/(4*mp.sqrt(mp.pi)*mp.gamma(1-s*(1-lam)))
        else:
            A=[half+s*lam,2*s*(1+lam),s*(2+lam),2*s,s*(1+lam)]
            B=[1+s*lam,1+2*s*lam,half+s*(2+lam),1+s*(1+lam)]
            pref=mp.sqrt(mp.pi)*s**(3-4*s)/(mp.mpf(4)**(1+s*lam)*(1+lam))
            pref*=mp.gamma(2*s*(1+lam))*mp.gamma(s*(2+lam))*mp.gamma(2*s)/(mp.gamma(1+s*lam)**2*mp.gamma(half+s*(2+lam)))
            x=pref*hyper_tail(A,B,*SETTINGS[working])
        return mp.nstr(x,working)

@lru_cache(None)
def compute(N,n,working):
    with mp.workdps(working):
        if N==1:
            assert n==2
            k=mp.root(3,3)*mp.gamma(mp.mpf(2)/3)/mp.gamma(mp.mpf(1)/3)
            x=k*k+1/k
        elif N==2:
            assert n==2
            x=mp.pi**2/8
        else:
            x=mp.mpf(boundary(N,n,1,working))+mp.mpf(boundary(N,n,-1,working))
        return mp.nstr(x,working)

class FullSpectralZeta(numberdb.Generator):
    table='T541'
    parameters=('N','n')
    type='R'
    digits=100
    rigour='heuristic (agreement-checked)'
    files=('generate.py','accelerate.py','checks.py','check.py','argument-selection.md','numerical-evidence.md','values.json','numerical-checks.json','skill-manifest.json')

    def enumerate(self):
        for N in range(1,13):
            for n in ((2,) if N<3 else (1,2)):
                yield {'N':N,'n':n}

    def value(self,params,digits):
        N,n=int(params['N']),int(params['n'])
        if N<1 or n not in (1,2) or (n==1 and N<3):
            raise ValueError('This generator implements the convergent first two moments only.')
        return numberdb.agreeing(lambda working:compute(N,n,working),at=(150,200))

if __name__=='__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
    gen=FullSpectralZeta()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1':
        print(gen.publish(message='Fill 22 convergent full oscillator traces with independent checks'))
    else:
        report=gen.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)