"""Numerical checks for irregular Coulomb values -- numberdb.org/T511. sage -python checks.py # all local checks, no credentials sage -python checks.py --stored # check database values; draft needs key Uses independent mpmath Coulomb G at 150 and 200 decimal digits; analytic U derivatives; M-based regular solutions with a product normalization; Bessel Y and cosine; and transcribed A&S / mpmath numerical controls. The reported JSON is reproducible evidence, not a certificate from mpmath. The value certification is from Arb enclosures and exact decimal endpoints. """ import numberdb.sage as numberdb from sage.rings.complex_arb import ComplexBallField from sage.rings.real_arb import RealBallField from sage.rings.real_mpfi import RealIntervalField from sage.rings.rational_field import QQ from sage.version import version as sage_version from numberdb._write import to_text from numberdb._generate import _stored_entries from generate import IrregularCoulomb from decimal import Decimal from fractions import Fraction from math import factorial import mpmath as mp import json import os import sys import time BITS = numberdb.bits(180) C = ComplexBallField(BITS) R = RealBallField(BITS) RI = RealIntervalField(BITS) def exact_decimal(text): return QQ(str(Fraction(Decimal(str(text).rstrip('?'))))) def interval_of(text): d = Decimal(str(text).rstrip('?')) center = QQ(str(Fraction(d))) ulp = QQ(10)**int(d.as_tuple().exponent) return RI(center-ulp, center+ulp), center, ulp def bounded(x, scale=1, places=130): assert x.is_finite(), 'nonfinite enclosure' assert x.rad_as_ball() < R(10)**(-places)*max(R(1),abs(R(scale))), 'wide enclosure' def zero(x, scale=1, places=130): bounded(x,scale,places) assert x.contains_zero(), 'identity residual excludes zero: %s' % x def regular(ell, eta, rho): """Independent regular solution and derivatives, using a product for C_l.""" e,r,i=C(QQ(eta)),C(QQ(rho)),C(0,1) if QQ(eta)==0: norm=C(1) else: norm=(2*C.pi()*e/((2*C.pi()*e).exp()-1)).sqrt() for j in range(1,ell+1): norm *= (j*j+e*e).sqrt() norm *= C(2)**ell / factorial(2*ell+1) a,b,z=ell+1-i*e,C(2*ell+2),2*i*r m=z.hypergeometric([a],[b]) m1=z.hypergeometric([a+1],[b+1]) m2=z.hypergeometric([a+2],[b+2]) pref=norm*r**(ell+1)*(-i*r).exp() A=(ell+1)/r-i f=pref*m df=pref*(A*m+2*i*a/b*m1) ddf=pref*((A*A-(ell+1)/(r*r))*m+4*i*A*a/b*m1-4*a*(a+1)/(b*(b+1))*m2) for x in (f,df,ddf): zero(x.imag(),abs(x.real())) return f.real(),df.real(),ddf.real() def irregular_derivatives(ell,eta,rho): """U derivatives and ODE residual independent of order recurrences.""" e,r,i=C(QQ(eta)),C(QQ(rho)),C(0,1) a,b,z=ell+1+i*e,C(2*ell+2),-2*i*r gamma=a.gamma() pref=(-i)**ell*(C.pi()*e/2).exp()*gamma/abs(gamma)*(i*r).exp()*z**(ell+1) u=z.hypergeometric_U(a,b) u1=z.hypergeometric_U(a+1,b+1) u2=z.hypergeometric_U(a+2,b+2) A=(ell+1)/r+i h=pref*u theta=r-e*(2*r).log()-ell*C.pi()/2+a.log_gamma().imag() h_dlmf=(i*theta).exp()*z**a*u zero((h-h_dlmf).real(),abs(h.real())) zero((h-h_dlmf).imag(),abs(h.imag())) dh=pref*(A*u+2*i*a*u1) ddh=pref*((A*A-(ell+1)/(r*r))*u+4*i*A*a*u1-4*a*(a+1)*u2) return h.real(),dh.real(),ddh.real(),h.imag(),dh.imag() def run(stored=False): started=time.monotonic() gen=IrregularCoulomb() if stored: doc=numberdb.table(gen.table) texts={tuple(identity.split(',')):text for identity,text in _stored_entries(doc).items()} expected={tuple(p[k] for k in gen.parameters) for p in gen.enumerate()} assert set(texts)==expected else: texts={} # Known nonzero controls precede comparisons; wrong sign must fail. g0=gen.value({'ell':'0','eta':'0','rho':'1'},100) zero(g0-R(1).cos()) assert not (g0+R(1).cos()).contains_zero() f0,df0,ddf0=regular(0,'0','1') zero(f0-R(1).sin());zero(df0-R(1).cos());zero(ddf0+R(1).sin()) _,dg0,_,_,_=irregular_derivatives(0,'0','1') w0=g0*df0-f0*dg0 zero(w0-1); assert not w0.contains_zero() # Controls for the decimal parser: 1.25 denotes [1.24,1.26]. control,center,ulp=interval_of('1.25') assert center==QQ(5)/4 and ulp==QQ(1)/100 assert QQ(124)/100 in control and QQ(123)/100 not in control counts={'decimal_enclosures':0,'mpmath_150':0,'mpmath_200':0,'precision_agreement':0, 'ode_G':0,'ode_F':0,'wronskian':0,'order_G':0,'order_F':0, 'derivative_up_G':0,'derivative_down_G':0,'bessel_Y':0,'cosine':0, 'handbook_values':0,'handbook_derivatives':0,'mpmath_printed':0, 'Hplus_imag_F':0} values={}; minbits=100000;maxrad=R(0); max_ulp_fraction=mp.mpf(0) handbook_g=('0.79445','1.0815','1.4969','2.0487') handbook_dg=('-0.67049','-0.60286','-0.56619','-0.79597') printed_g={1:'1.08148276',2:'1.496877075',3:'2.048694714'} printed_mdg={1:'0.6028279961',2:'0.5661803178',3:'0.7959909551'} for index,p in enumerate(gen.enumerate(),1): key=tuple(p[k] for k in gen.parameters) ell=int(p['ell']);e=R(QQ(p['eta']));r=R(QQ(p['rho'])) local=gen.value(p,100) bounded(local,local) minbits=min(minbits,int(local.accuracy()));maxrad=max(maxrad,local.rad_as_ball()) written=texts[key] if stored else to_text(local,digits=100) enclosure,center,ulp=interval_of(written) local_interval=RI(local) assert local_interval.lower()>=enclosure.lower() and local_interval.upper()<=enclosure.upper() counts['decimal_enclosures']+=1 g=R(enclosure) if stored else local scale=max(R(1),abs(g)) places=90 if stored else 125 values[key]=g f,df,ddf=regular(ell,p['eta'],p['rho']) gu,dg,ddg,imag_h,imag_dh=irregular_derivatives(ell,p['eta'],p['rho']) zero(g-gu,scale,places) zero(imag_h-f,max(abs(f),scale),125);zero(imag_dh-df,max(abs(df),scale),125) counts['Hplus_imag_F']+=1 coeff=1-2*e/r-ell*(ell+1)/(r*r) zero(ddg+coeff*g,scale*max(1,abs(coeff)),places);counts['ode_G']+=1 zero(ddf+coeff*f,max(1,abs(f)*abs(coeff)),125);counts['ode_F']+=1 zero(g*df-f*dg-1,max(1,abs(g*df),abs(f*dg)),places);counts['wronskian']+=1 def rr(j): return (1+e*e/(j*j)).sqrt() def ss(j): return j/r+e/j pup={'ell':str(ell+1),'eta':p['eta'],'rho':p['rho']} upkey=tuple(pup[k] for k in gen.parameters) if stored and upkey in texts: gp=R(interval_of(texts[upkey])[0]) else: gp=gen.value(pup,100) fp=regular(ell+1,p['eta'],p['rho'])[0] zero(dg-ss(ell+1)*g+rr(ell+1)*gp,max(scale*abs(ss(ell+1)),abs(gp)),places) counts['derivative_up_G']+=1 if ell: downkey=(str(ell-1),p['eta'],p['rho']) gm=values[downkey] fm=regular(ell-1,p['eta'],p['rho'])[0] zero(rr(ell)*gm-(ss(ell)+ss(ell+1))*g+rr(ell+1)*gp,max(abs(gm),abs(gp),scale*abs(ss(ell)+ss(ell+1))),places) counts['order_G']+=1 zero(rr(ell)*fm-(ss(ell)+ss(ell+1))*f+rr(ell+1)*fp,max(1,abs(fm),abs(fp)),125) counts['order_F']+=1 zero(dg-rr(ell)*gm+ss(ell)*g,max(scale*abs(ss(ell)),abs(gm)),places) counts['derivative_down_G']+=1 if p['eta']=='0': b=-(R.pi()*r/2).sqrt()*C(r).bessel_Y(C(QQ(ell)+QQ(1)/2)).real() bounded(b,scale) zero(g-b,scale,places);counts['bessel_Y']+=1 if ell==0:zero(g-r.cos(),scale,places);counts['cosine']+=1 if p['eta']=='2' and p['rho']=='5': assert abs(gu-R(exact_decimal(handbook_g[ell])))