"""Real zeros of the derivative of the Scorer function Gi -- numberdb.org/T517. Run with SageMath (or passagemath with numberdb): sage -pip install numberdb sage -python generate.py # verify every stored entry sage -python generate.py --fill # fill the existing draft Set NUMBERDB_API_KEY for authenticated draft access. --fill sends entries; it does not publish the table or offer it for review. The quantity is the signed location g'_n, not the height of Gi. Index zero selects the unique positive real zero; 1,...,12 select the negative zeros nearest zero, in decreasing order. DLMF 9.12 fixes the normalization. Certification uses 700-bit Arb arithmetic, 900 Maclaurin coefficients, and geometric tails on |x| <= 16 for Gi and its first two derivatives. For a term indexed by k, the ratios three indices later for derivative orders 0,1,2 are respectively R^3/((k+3)(k+2)), R^3/(k(k+2)), and R^3/(k(k-1)). After the last term in each residue class the ratios decrease. The absolute omitted tail is at most the last absolute term times q/(1-q). To count roots, adjacent rational cells cover the path from zero to the 12th negative root and [0,1]. On a cell of half-width h and midpoint m, Gronwall applied to (Gi,Gi') gives the uniform bound B = (|Gi(m)| + |Gi'(m)| + h/pi) exp(16 h) on both components, since max(1,|x|) <= 16. Thus Gi(x) is in Gi(m) +/- h B, and Gi'(x) is in Gi'(m) +/- h (R B + 1/pi), where R bounds |x| on the cell. Gi''(x) = x Gi(x) - 1/pi then certifies monotonicity of Gi'. Each cell either excludes a zero of Gi', or has a fixed sign of Gi'' and endpoint signs that determine whether it contains exactly one zero. Unresolved cells are subdivided. This certifies the negative-root indices; DLMF 9.12(ix) establishes uniqueness on the entire nonnegative axis. mpmath supplies candidates only. Every root is certified independently of the candidate accuracy by opposite ball signs at exact rational endpoints 2*10^-120 apart and fixed-sign Gi'' throughout that bracket. Root balls have about 120 decimal digits, leaving a measured 20-digit margin over the 100 significant digits stored. checks.py supplies independent integral and mpmath comparisons, with controls, and can repeat them on stored values. """ import os import sys from functools import lru_cache import numberdb.sage as numberdb from sage.rings.rational_field import QQ from sage.rings.real_arb import RealBallField import mpmath as mp PRECISION_BITS = 700 TERMS = 900 DOMAIN_RADIUS = 16 BRACKET_PLACES = 120 R = RealBallField(PRECISION_BITS) def sign(x): assert x.is_finite() return 1 if x > 0 else -1 if x < 0 else 0 def widen(x, error): return x + R(0, R(error).upper()) @lru_cache(None) def coefficients(): c = [1 / (R(3)**(QQ(7)/6) * R(QQ(2)/3).gamma()), 1 / (R(3)**(QQ(5)/6) * R(QQ(1)/3).gamma()), -1 / (2*R.pi())] for k in range(TERMS-3): c.append(c[k] / ((k+3)*(k+2))) return c @lru_cache(None) def derivative_coefficients(d): c = coefficients() if d == 0: return c if d == 1: return [k*c[k] for k in range(1, TERMS)] if d == 2: return [k*(k-1)*c[k] for k in range(2, TERMS)] raise ValueError('only derivative orders zero through two are implemented') @lru_cache(None) def tail(d): c = coefficients() total = R(0) for k in range(TERMS-3, TERMS): factor = [1, k, k*(k-1)][d] denominator = [(k+3)*(k+2), k*(k+2), k*(k-1)][d] q = R(DOMAIN_RADIUS)**3 / denominator assert 0 < q < 1 total += abs(c[k])*factor*R(DOMAIN_RADIUS)**(k-d)*q/(1-q) assert total < R('1e-400') return total def evaluate(x, d=0): x = R(x) assert abs(x) <= DOMAIN_RADIUS answer = R(0) for c in reversed(derivative_coefficients(d)): answer = answer*x+c answer = widen(answer, tail(d)) assert answer.is_finite() return answer @lru_cache(None) def point(x): return evaluate(x), evaluate(x, 1) def local_enclosures(a, b): m, h = (a+b)/2, (b-a)/2 y, v = point(m) bound = (abs(y)+abs(v)+R(h)/R.pi())*(16*R(h)).exp() y_box = widen(y, R(h)*bound) v_box = widen(v, R(h)*(max(abs(a),abs(b))*bound+1/R.pi())) acceleration = R(m, h)*y_box-1/R.pi() return v_box, acceleration def classify(a, b, records, depth=0): v_box, acceleration = local_enclosures(a,b) if sign(v_box): records.append({'a':str(a),'b':str(b),'test':'excludes-zero'}) return [] if sign(acceleration): sa, sb = sign(point(a)[1]), sign(point(b)[1]) assert sa and sb roots = sa != sb records.append({'a':str(a),'b':str(b),'test':'monotone', 'left_sign':sa,'right_sign':sb,'root':roots}) return [(a,b)] if roots else [] assert depth < 20, 'interval root count did not resolve' m = (a+b)/2 return classify(a,m,records,depth+1)+classify(m,b,records,depth+1) @lru_cache(None) def root_cells(): records=[] positive=[] for k in range(64): positive += classify(QQ(k)/64,QQ(k+1)/64,records) assert len(positive)==1 negative=[] for k in range(1024): negative += list(reversed(classify(-QQ(k+1)/64,-QQ(k)/64,records))) if len(negative)==12: break assert len(negative)==12 return positive+negative, records def mp_rational(q): return mp.mpf(int(q.numerator()))/int(q.denominator()) @lru_cache(None) def certified_root(n): cells, _ = root_cells() a,b = cells[n] with mp.workdps(170): z = mp.findroot(lambda x:mp.diff(mp.scorergi,x), (mp_rational(a),mp_rational(b))) center = QQ(int(mp.nint(z*mp.mpf(10)**BRACKET_PLACES)))/10**BRACKET_PLACES epsilon = QQ(1)/10**BRACKET_PLACES lo,hi = center-epsilon,center+epsilon assert a < lo < hi < b assert sign(evaluate(lo,1))*sign(evaluate(hi,1)) == -1 assert sign(local_enclosures(lo,hi)[1]) ball=R(center,epsilon) assert ball.is_finite() and R(ball.rad()) < R('1.01e-120') return ball,lo,hi class ScorerGiPrimeZeros(numberdb.Generator): table='T517' parameters=('n',) type='R' digits=100 rigour='proven' files=('generate.py','checks.py','check.py') def enumerate(self): for n in range(13): yield {'n':str(n)} def value(self, params, digits): if digits > 100: raise ValueError('increase and revalidate the certificate for more than 100 digits') n=int(params['n']) if not 0 <= n <= 12: raise ValueError('the certified index range is zero through twelve') value=certified_root(n)[0] if n==0: return {'number':value,'comment':'The unique positive real zero of $\\operatorname{Gi}^{\\prime}$. It is the location of a local maximum of $\\operatorname{Gi}$.'} return value if __name__=='__main__': if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1': os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip() generator=ScorerGiPrimeZeros() if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1': print(generator.publish(message='Fill thirteen certified signed Gi derivative zeros')) else: report=generator.verify(sample=None) print(report) sys.exit(0 if report.ok else 1)