"""RIGOROUS certificate for hypothesis (ii) of Polymath's Theorem 1.2 (arXiv:1904.12438, Thm 1.2/1.3, (20)-(24)): H_{t0}(x+iy) != 0 for x >= X + sqrt(1-y0^2), y0 <= y <= sqrt(1-2 t0). Method (see the accompanying paper): optimised mollifier M = sum_{d<=D} lam_d d^{-s*}, lam_1 = 1 (lam_d are fixed float64 numbers, hence exact rationals). For every point, |M| (|A| - |gamma||B|) = |MA| - |gamma||M Bbar| >= LB, |f_t0| >= LB / sum|lam_d| d^{-Re s*} > eA+eB+eC0. Transcendentals: Arb balls (python-flint). Bulk sums: float64 with outward rounding (see common.py). Usage: python cert_region_ii.py lam.npy 0.15499 0.2 5998313329980.3135 """ import sys, time, math import numpy as np import os as _os; sys.path.insert(0,_os.path.dirname(_os.path.abspath(__file__))) from common import arb, init, lo, hi, dn, up, sum_up, A, geom_blocks, table_b, U import os HERE=os.path.dirname(os.path.abspath(__file__)); sys.path.insert(0,HERE) if len(sys.argv)<5: sys.exit('usage: python cert_region_ii.py WEIGHTS.npy t0 y0 X [ystart h] (normally called by certify_all.py)') T0s,Y0s,Xs=sys.argv[2],sys.argv[3],sys.argv[4] lam=np.load(sys.argv[1]).astype(np.float64); Dn=len(lam); D=np.arange(1,Dn+1) assert lam[0]==1.0 PI=init(96) T0=A(T0s); Y0=A(Y0s); X=A(Xs) t=T0; Xii=X+(1-Y0*Y0).sqrt(); YTOP=(1-2*T0).sqrt() t_f=float(T0s) LOGD=[arb(int(d)).log() for d in D]; ALAM=[arb(float(abs(l))) for l in lam] def sig_lb(x,y): # Polymath (21); the bracket (...)_+ is <= 2 for x>=200, 0<=y<=1 return (1+y)/2+(t/4)*(x/(4*PI)).log()-t/(x*x) def kap_ub(x,y): return t*y/(2*(x-6)) # (22) def gam_ub(x,y): return (arb('0.02')*y).exp()*((x/(4*PI)).log()*(-y/2)).exp() # (20) def Nof(x): s=(x/(4*PI)).sqrt()+t/16; a,b=math.floor(lo(s)),math.floor(hi(s)) assert a==b, 'N(x) ambiguous'; return a def xN(m): return 4*PI*(arb(m)-t/16)**2 # ---------------------------------------------------------------- b_m table (filled lazily) BT={'N':0} def need_b(N): if N>BT['N']: t0=time.time(); N2=max(N,int(BT['N']*1.25)); blo,bhi=table_b(N2,t,BT['N']+1) if BT['N']==0: BT.update(N=N2,blo=blo,bhi=bhi) else: BT.update(N=N2,blo=np.concatenate([BT['blo'],blo]),bhi=np.concatenate([BT['bhi'],bhi])) N=N2 print(' [b_m table up to %d: %.0fs]'%(N,time.time()-t0),flush=True) def rmul_pos(a_lo,a_hi,b_lo,b_hi): # interval product, b>=0 p=np.stack([a_lo*b_lo,a_lo*b_hi,a_hi*b_lo,a_hi*b_hi]); return dn(p.min(0)),up(p.max(0)) def main_sum(xa,Na,ya,yb,chunk=4000000): need_b(Na); blo=BT['blo'][:Na]; bhi=BT['bhi'][:Na] sA=lo(sig_lb(xa,ya)) eBa=ya-sig_lb(xa,ya)+kap_ub(xa,ya); eBb=yb-sig_lb(xa,yb)+kap_ub(xa,yb) eB=max(hi(eBa),hi(eBb)); assert eB<0 g=hi(gam_ub(xa,ya)) dyl=np.array([lo((-ya*LOGD[j]).exp()) for j in range(Dn)]); dyh=np.array([hi((-ya*LOGD[j]).exp()) for j in range(Dn)]) ddy=np.array([hi((-ya*LOGD[j]).exp()-(-yb*LOGD[j]).exp()) for j in range(Dn)]) cl=np.where(lam>=0,lam*dyl,lam*dyh); ch=np.where(lam>=0,lam*dyh,lam*dyl) # lam_j d^{-ya} interval cl=dn(cl); ch=up(ch) nmax=Dn*Na; Bk=geom_blocks(nmax); nb=len(Bk)-1 sA_arb=arb(sA); eB_arb=arb(eB) wA=np.array([hi((-sA_arb*arb(int(b)).log()).exp()) for b in Bk[:-1]]) # n^{-sA} <= B_k^{-sA} wB=np.array([hi((eB_arb*arb(int(b)).log()).exp()) for b in Bk[:-1]]) # n^{eB} <= B_k^{eB} (eB<0) accA=np.zeros(nb); accB=np.zeros(nb) for lo_n in range(1,nmax+1,chunk): hi_n=min(lo_n+chunk,nmax+1); L=hi_n-lo_n Al=np.zeros(L); Ah=np.zeros(L); Bl=np.zeros(L); Bh=np.zeros(L); E=np.zeros(L) for j in range(Dn): d=j+1; m0=max(1,-(-lo_n//d)); m1=min(Na,(hi_n-1)//d) if m1=0: vl,vh=dn(lj*bl),up(lj*bh) else: vl,vh=dn(lj*bh),up(lj*bl) Al[r]=dn(Al[r]+vl); Ah[r]=up(Ah[r]+vh) ul,uh=rmul_pos(cl[j],ch[j],bl,bh); Bl[r]=dn(Bl[r]+ul); Bh[r]=up(Bh[r]+uh) E[r]=up(E[r]+up(up(abs(lj)*up(ddy[j]+up(1e-12*dyh[j])))*bh)) # + d^{conj kappa}-1 <= 1e-12 (|kappa|<1e-13, d<=240) absA=np.maximum(np.abs(Al),np.abs(Ah)); absB=up(np.maximum(np.abs(Bl),np.abs(Bh))+E) if lo_n==1: absA[0]=0.0; assert Al[0]<=1.0<=Ah[0] # c_1 = lam_1 b_1 = 1 is the leading term n=np.arange(lo_n,hi_n); k=np.searchsorted(Bk,n,side='right')-1 gm=1+2*L*U accA+=up(np.bincount(k,weights=absA,minlength=nb)*gm); accB+=up(np.bincount(k,weights=absB,minlength=nb)*gm) accA=up(accA); accB=up(accB) SA=sum_up(up(accA*wA)); SB=sum_up(up(accB*wB)) LB=float(dn(dn(1-SA)-up(g*SB))) Mmax=sum_up(np.array([hi(ALAM[j]*(-sA_arb*LOGD[j]).exp()) for j in range(Dn)])) return LB,Na,Mmax def F_arb(v,ya,yb): """upper bound for all new terms with index v (v real >= 1): A-part + B-part, evaluated at x_v (decreasing in v)""" xv=xN(v); Lv=arb(v).log() if isinstance(v,int) else v.log() s=sig_lb(xv,ya); eB=ya-sig_lb(xv,ya)+kap_ub(xv,ya); eB2=yb-sig_lb(xv,yb)+kap_ub(xv,yb) eBu=arb(max(hi(eB),hi(eB2))) bm=(t/4*Lv*Lv).exp(); g=gam_ub(xv,ya) SA=arb(0); SB=arb(0) for j in range(Dn): if lam[j]==0: continue SA+=ALAM[j]*(-s*LOGD[j]).exp(); SB+=ALAM[j]*((eBu-ya)*LOGD[j]).exp() SB=SB*(1+arb('1e-12')) # d^{conj kappa} factor: |d^{conj kappa}| <= D^{|kappa|} <= 1+1e-12 return hi(bm*(-s*Lv).exp()*SA + g*bm*(eBu*Lv).exp()*SB) def extra(Na,Nb,ya,yb,r=1+2.0**-9): if Nb<=Na: return 0.0 starts=[Na+1]; v=float(Na+1) while True: v=max(v*r,starts[-1]+1); s=int(math.ceil(v)) if s>Nb: break starts.append(s) starts.append(Nb+1); tot=[] for a,b in zip(starts[:-1],starts[1:]): tot.append(up((b-a)*F_arb(a,ya,yb))) return sum_up(np.array(tot)) def future_tail(M,ya,yb,r=1.01,LMAX=60): """sum_{m>M} F(m) <= int_M^inf F <= sum_k (v_{k+1}-v_k) F(v_k) + e^{-LMAX}; for log v >= LMAX one has F(v) <= v^{-2} (checked at v=e^{LMAX}; d/dL[log F + 2L] <= -(t/2)L + 2 < 0 there)""" v=arb(M); tot=[] while hi(v.log())2.5 return up(sum_up(np.array(tot))+hi(1/VE)) def float_extra_est(Na,Nb,ya): # cheap planning estimate only (not part of the certificate) K=min(Nb-Na,20000); m=np.linspace(Na+1,Nb,K); wgt=(Nb-Na)/K; x=4*math.pi*(m-t_f/16)**2; L=np.log(m) s=(1+ya)/2+(t_f/4)*np.log(x/(4*math.pi)); b=np.exp(t_f/4*L*L) return float(wgt*np.sum(b*np.exp(-s*L))*np.sum(np.abs(lam)*D**(-s[0]))*1.25) def errors_majorant(ya_all,yb_all): """uniform bound for eA+eB+eC0 over x>=Xii, y in [ya_all,yb_all] (all pieces non-increasing in x)""" x=Xii; N=Nof(x); lx=(x/(4*PI)).log() eC=((-(1+ya_all)/4)*lx).exp()*(-(t/16)*lx*lx+arb('1.24')*(3*yb_all+1)/(N-arb('0.125'))+(3*(lx*lx+PI*PI/4).sqrt()+arb('10.44'))/(x-12)).exp() # sum_{n<=N(x)} b_n n^{-Re s*} for every x>=Xii <= (value at Xii) + (all later terms) need_b(N); s=lo(sig_lb(x,ya_all)); Bk=geom_blocks(N) w=np.array([hi((-arb(s)*arb(int(b)).log()).exp()) for b in Bk[:-1]]) acc=np.zeros(len(Bk)-1); n=np.arange(1,N+1); k=np.searchsorted(Bk,n,side='right')-1 acc=up(np.bincount(k,weights=BT['bhi'][:N],minlength=len(acc))*(1+2*N*U)) S=sum_up(up(acc*w)) lam_save=None fac=hi(((t*t/16)*lx*lx+arb('0.626'))/(x-arb('6.66'))) fac=hi((arb(fac)).exp()-1) S=up(S+future_tail(N,ya_all,yb_all)) # later terms: F_arb >= plain terms since lam_1=1 eAB=hi((1+arb('0.02').exp()*arb('1.0001'))*arb(S)*arb(fac)) return hi(eAB+eC) if __name__=='__main__': T=time.time() print('Certificate (ii): t0=%s y0=%s X=%s D=%d'%(T0s,Y0s,Xs,Dn),flush=True) ytop=hi(YTOP); E=errors_majorant(Y0,arb(ytop)) # ytop: float upper bound of sqrt(1-2 t0) print('eA+eB+eC0 <= %.3e for all x>=X+sqrt(1-y0^2), y in [y0,%.6f]'%(E,ytop),flush=True) print('YTOP_HI %r (float upper bound of sqrt(1-2 t0); the y-intervals must reach it)'%ytop,flush=True) ya=Y0; ya_f=float(Y0s); h=0.005; worst=1e9; nint=0 if len(sys.argv)>6: # resume at y=ystart (must be <= the last certified y) with initial width h ya_f=float(sys.argv[5]); h=float(sys.argv[6]); ya=arb(ya_f); assert ya_f>=float(Y0s) print('RESUMING at y=%.6f (intervals below were certified in an earlier run)'%ya_f,flush=True) while ya_f=0.9*LBa and Nb>Na: Nb=Na+(Nb-Na)//2; ex=extra(Na,Nb,ya,yb) assert ex halving'%(ya_f,yb_f),flush=True); h/=2; assert h>1e-5; continue lbs=[float(dn(dn(r[2]-r[3])/r[4])) for r in rows]; mlb=min(lbs); assert mlb>E worst=min(worst,mlb); nint+=1 print('OK y in [%r,%r]: %d chunks, certified |f_t0| >= %.5f [%.0fs]'%(ya_f,yb_f,len(rows),mlb,time.time()-T),flush=True) for r in rows: print(' x_a~%.6e N_a=%d LB>=%.5f extra/tail<=%.5f |M|<=%.4f'%r,flush=True) ya=yb; ya_f=yb_f; h*=2 yfrom=float(sys.argv[5]) if len(sys.argv)>6 else float(Y0s) print('CERTIFIED (ii) [Arb]: t0=%s y0=%s X=%s : |f_t0| >= %.5f > eA+eB+eC0 (<= %.2e) for all x>=X+sqrt(1-y0^2), y in [%.6f, %.6f]; %d y-intervals; %.0fs'%(T0s,Y0s,Xs,worst,E,yfrom,ytop,nint,time.time()-T))