"""Real zeros of the Scorer function Gi -- numberdb.org/T516. Reproduce with SageMath and numberdb (mpmath is included in SageMath): sage -pip install numberdb sage -python generate.py # compare all twelve stored entries sage -python generate.py --check # local certification, no network sage -python generate.py --fill # fill this draft, NUMBERDB_API_KEY set This does not publish or offer a draft for review. Every computation establishes both the enclosures and the indexing anew. mpmath only proposes narrow brackets; Arb certifies them. The first twelve negative roots are enumerated from zero. At 700 bits, a 900-term Maclaurin polynomial encloses Gi and its derivative on [-16,0]. Each of the three residue-class tails is bounded geometrically, using its first omitted term and its decreasing ratio. For derivative order j=0,1,2, the ratios are respectively R^3/((k+3)(k+2)), R^3/((k+2)k), R^3/(k(k-1)). For indexing, on [m-h,m+h], Gronwall applied to (y,y') with the maximum norm gives M=(|y(m)|+|y'(m)|+h/pi)*exp(16h) as a bound for both |y| and |y'|. Thus y(m)+[-hM,hM] and y'(m)+[-h(16M+1/pi),h(16M+1/pi)] enclose y and y'. Each segment either excludes zero in y, or has derivative of fixed sign and its endpoint signs determine whether it has zero or one root. Inconclusive segments are bisected. No homogeneous Sturm theorem is assumed. """ 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 BITS = 700 TERMS = 900 RADIUS = 16 ROOT_RADIUS = QQ(1) / 10**120 RF = RealBallField(BITS) def rational_decimal(s): """Parse a finite decimal exactly, without a binary floating conversion.""" s = str(s).lower() mantissa, _, exponent = s.partition('e') exponent = int(exponent or 0) whole, dot, frac = mantissa.partition('.') return QQ(int(whole + frac)) * QQ(10)**(exponent-len(frac)) @lru_cache(None) def coefficients(): c = [RF(1)/(RF(3)**(QQ(7)/6)*RF(QQ(2)/3).gamma()), RF(1)/(RF(3)**(QQ(5)/6)*RF(QQ(1)/3).gamma()), -RF(1)/(2*RF.pi())] for k in range(TERMS): c.append(c[k]/((k+3)*(k+2))) return c @lru_cache(None) def polynomial_and_tail(order): c = coefficients() def factor(k): f = 1 for j in range(order): f *= k-j return f p = [c[k]*factor(k) for k in range(order, TERMS)] tail = RF(0) for k in range(TERMS, TERMS+3): if order == 0: denominator = (k+3)*(k+2) elif order == 1: denominator = (k+2)*k elif order == 2: denominator = k*(k-1) else: raise ValueError('supported derivative orders: 0, 1, 2') q = RF(RADIUS)**3/denominator assert q < 1 tail += abs(c[k])*factor(k)*RF(RADIUS)**(k-order)/(1-q) assert tail.is_finite() and tail < RF(10)**-160 return p, tail def gi(x, order=0): x = RF(x) assert abs(x) <= RADIUS p, tail = polynomial_and_tail(order) result = RF(0) for c in reversed(p): result = result*x+c result = result.add_error(tail.upper()) assert result.is_finite() return result def sign(y): if y > 0: return 1 if y < 0: return -1 return 0 def local_enclosures(a, b): m, h = (a+b)/2, (b-a)/2 y, yp = gi(m), gi(m, 1) M = (abs(y)+abs(yp)+RF(h)/RF.pi())*(RF(RADIUS)*h).exp() ybox = y.add_error((RF(h)*M).upper()) dbox = yp.add_error((RF(h)*(RADIUS*M+1/RF.pi())).upper()) assert ybox.is_finite() and dbox.is_finite() return ybox, dbox @lru_cache(None) def certify(): """Return twelve rational brackets and a full finite-interval root count.""" stats = {'excluded_by_value': 0, 'excluded_by_monotonicity': 0, 'root_segments': 0, 'subdivisions': 0} coarse = [] def visit(a, b, depth=0): assert depth < 30 ybox, dbox = local_enclosures(a,b) if sign(ybox): stats['excluded_by_value'] += 1 return if sign(dbox): sa, sb = sign(gi(a)), sign(gi(b)) assert sa and sb if sa != sb: coarse.append((a,b)) stats['root_segments'] += 1 else: stats['excluded_by_monotonicity'] += 1 return stats['subdivisions'] += 1 m = (a+b)/2 visit(m,b,depth+1) visit(a,m,depth+1) for i in range(1, 16*16+1): a,b = -QQ(i)/16,-QQ(i-1)/16 visit(a,b) if len(coarse) >= 12: assert len(coarse) == 12 stats['left_endpoint'] = str(a) break assert len(coarse) == 12 brackets=[] for a,b in coarse: with mp.workdps(170): ma,mb = mp.mpf(str(a.numerator()))/int(a.denominator()),mp.mpf(str(b.numerator()))/int(b.denominator()) candidate=mp.findroot(mp.scorergi,(ma,mb)) mid=rational_decimal(mp.nstr(candidate,160)) lo,hi=mid-ROOT_RADIUS,mid+ROOT_RADIUS assert a < lo < hi < b assert sign(gi(lo))*sign(gi(hi)) == -1 assert sign(local_enclosures(lo,hi)[1]) brackets.append((lo,hi)) assert all(brackets[i][0] > brackets[i+1][1] for i in range(11)) return brackets,stats class ScorerGiZeros(numberdb.Generator): files = ('generate.py', 'checks.py') table = 'T516' parameters = ('n',) type = 'R' digits = 100 rigour = 'proven' def enumerate(self): for n in range(1,13): yield {'n':n} def value(self,params,digits): if digits > 100: raise ValueError('This certificate is configured for 100 digits.') n=int(params['n']) if not 1 <= n <= 12: raise ValueError('The certified prefix has twelve roots.') a,b=certify()[0][n-1] return RF((a+b)/2).add_error(RF((b-a)/2).upper()) if __name__ == '__main__': if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1': os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip() g=ScorerGiZeros() if '--check' in sys.argv: print(certify()[1]) for p in g.enumerate(): print(p,g.value(p,100)) elif '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1': print(g.publish(message='Certify the first twelve signed real zeros, including their indices')) else: report=g.verify(sample=None) print(report) sys.exit(0 if report.ok else 1)