"""RIGOROUS certificate for hypothesis (iii) of Polymath's Theorem 1.2 (barrier):
   H_t(x+iy) != 0 for X <= x <= X+sqrt(1-y0^2), sqrt(y0^2+2(t0-t)) <= y <= sqrt(1-2t), 0 <= t <= t0.
For each time step [t1,t2] we use the rectangle R = [X, X+1] x [ylo, 1], ylo = sqrt(y0^2+2(t0-t2)), which contains
the region for every t in the step.  With Dz, Dt (Polymath Lemma 8.4) and E (eqs. (23),(24)):
   min_j |f_t1(z_j)| - Dz*h/2 - Dt*(t2-t1) - E > 0  and  min_j |f_t1(z_j)| > Dz*h,  winding number 0
imply (Rouche) that H_t has no zero in R for all t in (t1,t2]  (N is constant on R, so f_t is holomorphic there).
t = 0 itself is covered by hypothesis (i) (RH verified to height 3e12 >= (X+1)/2).
f_t is evaluated by Polymath's Taylor scheme (Sec. 7) on dyadic blocks, with complex-interval arithmetic.
Usage: python cert_barrier.py 5998313329980.3135 0.15499 0.2
"""
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, acb, init, lo, hi, dn, up, sum_up, A, geom_blocks, U
if len(sys.argv)<4:
    sys.exit('usage: python cert_barrier.py X t0 y0   (normally called by certify_all.py)')
Xs,T0s,Y0s=sys.argv[1],sys.argv[2],sys.argv[3]
PI=init(128)
X=A(Xs); XC=X+arb('0.5'); XR=X+1; T0=A(T0s); Y0=A(Y0s); t0f=float(T0s); y0f=float(Y0s)
YC=arb('0.6'); ycf=0.6
s_=(X/(4*PI)).sqrt(); s2_=(XR/(4*PI)).sqrt()+T0/16
N=math.floor(lo(s_)); assert math.floor(hi(s2_))==N, 'N not constant on the barrier'
print('Barrier certificate: X=%s t0=%s y0=%s N=%d'%(Xs,T0s,Y0s,N),flush=True)
AMAX,BMAX=44,10
# ------------------------------------------------------------------ per-n data (Arb)
T=time.time()
blocks=[]; k=0
while 2**k<=N: blocks.append((2**k,min(2**(k+1)-1,N))); k+=1
TWO_PI=2*PI; HALF=XC/2
def interval_sum(lo_arr,hi_arr):
    n=lo_arr.size; g=n*U/(1-n*U)*2
    s_lo=float(np.sum(lo_arr)); s_hi=float(np.sum(hi_arr))
    return float(dn(dn(s_lo)-g*float(np.sum(np.abs(lo_arr)))*1.01)), float(up(up(s_hi)+g*float(np.sum(np.abs(hi_arr)))*1.01))
def imul(al,ah,bl,bh):
    p=np.stack([al*bl,al*bh,ah*bl,ah*bh]); return dn(p.min(0)),up(p.max(0))
MOM=[]
for (n0,n1) in blocks:
    m=n1-n0+1; cl=np.empty(m); ch=np.empty(m); sl=np.empty(m); sh=np.empty(m)
    a1l=np.empty(m); a1h=np.empty(m); a2l=np.empty(m); a2h=np.empty(m); ul=np.empty(m); uh=np.empty(m)
    l0=arb(n0).log()
    for i,n in enumerate(range(n0,n1+1)):
        L=arb(n).log(); v=HALF*L; q=math.floor(float((v/TWO_PI).mid())); th=v-TWO_PI*q
        c=th.cos(); s=th.sin()
        cl[i],ch[i]=lo(c),hi(c); sl[i],sh[i]=lo(s),hi(s)
        e1=((-(1+YC)/2)*L).exp(); e2=((-(1-YC)/2)*L).exp()
        a1l[i],a1h[i]=lo(e1),hi(e1); a2l[i],a2h[i]=lo(e2),hi(e2)
        u=L-l0; ul[i],uh[i]=max(0.0,lo(u)),hi(u)
    # A-moments: sum a1 e^{+i th} u^j ; B-moments: sum a2 e^{-i th} u^j
    Ar=imul(a1l,a1h,cl,ch); Ai=imul(a1l,a1h,sl,sh); Br=imul(a2l,a2h,cl,ch); Bi=imul(a2l,a2h,-sh,-sl)
    J=AMAX+2*BMAX; pl=np.ones(m); ph=np.ones(m)
    mom=np.zeros((J+1,8))
    for j in range(J+1):
        row=[]
        for (zl,zh) in (Ar,Ai,Br,Bi):
            vl,vh=imul(zl,zh,pl,ph); row+=list(interval_sum(vl,vh))
        mom[j]=row
        pl=dn(pl*ul); ph=up(ph*uh)
    absA=sum_up(a1h); absB=sum_up(a2h); umax=float(uh.max())
    MOM.append(dict(n0=n0,l=l0,mom=mom,absA=absA,absB=absB,umax=umax))
print('moments: %d blocks, %.0fs'%(len(MOM),time.time()-T),flush=True)
FA=[1/math.factorial(a) for a in range(AMAX+1)]; FB=[1/math.factorial(b) for b in range(BMAX+1)]
# ------------------------------------------------------------------ complex interval helpers (numpy rectangles)
def cmul(a,b):
    rr=imul(a[0],a[1],b[0],b[1]); ii=imul(a[2],a[3],b[2],b[3]); ri=imul(a[0],a[1],b[2],b[3]); ir=imul(a[2],a[3],b[0],b[1])
    return (dn(rr[0]-ii[1]),up(rr[1]-ii[0]),dn(ri[0]+ir[0]),up(ri[1]+ir[1]))
def cadd(a,b): return (dn(a[0]+b[0]),up(a[1]+b[1]),dn(a[2]+b[2]),up(a[3]+b[3]))
def from_acb(lst):
    return (np.array([lo(z.real) for z in lst]),np.array([hi(z.real) for z in lst]),np.array([lo(z.imag) for z in lst]),np.array([hi(z.imag) for z in lst]))
def coef(mom,a,tq,part):   # r_a = sum_b (t/4)^b/b! m_{a+2b} / a!  as complex interval (scalars), outward-rounded
    o=0 if part=='A' else 4
    re_l=re_h=im_l=im_h=0.0; ab=0.0
    for b in range(BMAX+1):
        w=tq[b]*FB[b]*FA[a]; row=mom[a+2*b]
        re_l+=w*row[o]; re_h+=w*row[o+1]; im_l+=w*row[o+2]; im_h+=w*row[o+3]
        ab+=w*(abs(row[o])+abs(row[o+1])+abs(row[o+2])+abs(row[o+3]))
    sl=8*(BMAX+4)*U*ab+1e-300     # rounding of w (factorials, powers) and of the short sums
    return (float(dn(re_l-sl)), float(up(re_h+sl)), float(dn(im_l-sl)), float(up(im_h+sl)))
def alpha(s): return 1/(2*s)+1/(s-1)+(s/(2*PI)).log()/2
def logM0(s): return s.log()+(s-1).log()-(s/2)*PI.log()+((2*PI).sqrt()/16).log()+(s/2-arb('0.5'))*(s/2).log()-s/2
def f_eval(dx,y,t):
    """complex-interval enclosure of f_t at x = XC+dx (dx, y exact floats); returns rectangle arrays"""
    tA=arb(t); npt=len(dx); tq=[1.0]
    for b_ in range(BMAX): tq.append(tq[-1]*(t/4))      # correctly rounded products (relative error covered by coef slack)
    assert float(tA)==t
    D_A=[];D_B=[];G=[]
    for i in range(npt):
        x=XC+arb(float(dx[i])); yy=arb(float(y[i]))
        s1=acb((1+yy)/2,-x/2); s2=acb((1-yy)/2,x/2)
        a1=alpha(s1); a2=alpha(s2)
        D_A.append(acb(arb(float(y[i]))-YC, -arb(float(dx[i])))/2+tA/2*a1)
        D_B.append(acb(YC-arb(float(y[i])), arb(float(dx[i])))/2+tA/2*a2)
        G.append(((tA/4)*(a2*a2-a1*a1)+logM0(s2)-logM0(s1)).exp())
    tot={}; trn={'A':np.zeros(npt),'B':np.zeros(npt)}
    fA=math.factorial(AMAX+1); fB=math.factorial(BMAX+1)
    for part,DL in (('A',D_A),('B',D_B)):
        acc_tot=(np.zeros(npt),np.zeros(npt),np.zeros(npt),np.zeros(npt))
        for B in MOM:
            l=B['l']
            pref=[((tA/4)*l*l-d*l).exp() for d in DL]; cc=[(tA/2)*l-d for d in DL]
            C=from_acb(cc); P=from_acb(pref)
            r=coef(B['mom'],AMAX,tq,part); acc=tuple(np.full(npt,v) for v in r)
            for a in range(AMAX-1,-1,-1):
                r=coef(B['mom'],a,tq,part); acc=cadd(cmul(acc,C),tuple(np.full(npt,v) for v in r))
            acc_tot=cadd(acc_tot,cmul(P,acc))
            # truncation of the a-series (a>AMAX) and of the b-series (b>BMAX)
            cmax=float(up(np.max(up(np.maximum(np.abs(C[0]),np.abs(C[1]))+np.maximum(np.abs(C[2]),np.abs(C[3]))))))
            cu=arb(cmax)*arb(B['umax']); w2=tA/4*arb(B['umax'])**2          # |c|u and (t/4)u^2, upper bounds (Arb)
            tr=hi((cu**(AMAX+1))/fA*(cu+w2).exp()+(w2**(BMAX+1))/fB*(w2+cu).exp())   # remainders of the a- and b-series
            pabs=up(np.maximum(np.abs(P[0]),np.abs(P[1]))+np.maximum(np.abs(P[2]),np.abs(P[3])))
            trn[part]=up(trn[part]+up(up(pabs*(B['absA'] if part=='A' else B['absB']))*tr))
        tot[part]=acc_tot
    Gi=from_acb(G)
    f=cadd(tot['A'],cmul(Gi,tot['B']))
    gabs=up(np.maximum(np.abs(Gi[0]),np.abs(Gi[1]))+np.maximum(np.abs(Gi[2]),np.abs(Gi[3])))
    err=up(up(trn['A']+up(gabs*trn['B']))+1e-300)
    return f,err
def absmin(f):
    rx=np.where((f[0]<=0)&(f[1]>=0),0.0,np.minimum(np.abs(f[0]),np.abs(f[1])))
    ry=np.where((f[2]<=0)&(f[3]>=0),0.0,np.minimum(np.abs(f[2]),np.abs(f[3])))
    return dn(np.sqrt(dn(rx*rx+ry*ry))*(1-1e-15))
# ------------------------------------------------------------------ bounds (Lemma 8.4, (23), (24)) via monotone blocks
BK=geom_blocks(N,r=1+2.0**-10)
BKl=[arb(int(b)).log() for b in BK[:-1]]; BKe=[arb(int(b)-1).log() if b>1 else arb(0) for b in BK[1:]]   # log of block start / block end
CNT=np.diff(BK).astype(float)
LX=(XR/(4*PI)).log(); LXLO=(X/(4*PI)).log()
def bounds(t1,t2,ylo):
    tA=arb(t1); tB=arb(t2); y=arb(ylo)
    rs=(1+y)/2+(tA/4)*LXLO-tA/(X*X)                    # (21) lower bound for Re s* on R at time t1
    g=(arb('0.02')).exp()*((-y/2)*LXLO).exp()*((tB/(2*(X-6)))*arb(N).log()).exp()   # |gamma| N^{|kappa|}, y<=1
    lam=(arb('0.5')+tB/(4*(X-6)))*(((4+XR*XR).sqrt()/(4*PI)).log()+PI+3/X)
    Dz=arb(0); Dt=arb(0); S=arb(0); SB=arb(0)
    for k in range(len(CNT)):
        Ls=BKl[k]; Le=BKe[k]; cnt=arb(CNT[k])
        w=((tA/4)*Ls*Ls-rs*Ls).exp()                    # b_n n^{-Re s*}: decreasing in n  -> value at block start
        wB=g*w*(y*Ls).exp()                              # g b_n n^{y-Re s*}: decreasing in n
        base=Le/4*(XR/(4*PI)).log()-Le*Ls/4+(PI/8)*Le+2*Le/(X-6)   # (1/4)L log(x/4pi n) + pi/8 L + ..., L<=Le, log n>=Ls
        Dz+=cnt*(w*(Le/2+tB*Le/(4*(X-6)))+wB*(tB*Le/(4*(X-6))+lam))
        Dt+=cnt*(w*base+wB*(base+arb('0.25')*(PI/2+8/(X-6))*(LX+8/(X-6))))
        S+=cnt*(w+wB)
    fac=(((tB*tB/16)*LX*LX+arb('0.626'))/(X-arb('6.66'))).exp()-1
    eC=((-(1+y)/4)*LXLO).exp()*(-(tA/16)*LXLO*LXLO+arb('1.24')*4/(N-arb('0.125'))+(3*(LX*LX+PI*PI/4).sqrt()+arb('10.44'))/(X-12)).exp()
    return hi(Dz),hi(Dt),hi(S*fac+eC)
def mesh(ylo,h):
    nx=int(math.ceil(1.0/h)); H=1.0-ylo; ny=int(math.ceil(H/h))
    xb=-0.5+np.arange(nx+1)/nx; yb=ylo+H*np.arange(ny+1)/ny; yb[-1]=1.0; yb[0]=ylo
    xs=np.concatenate([xb,np.full(ny,0.5),xb[::-1][1:],np.full(ny-1,-0.5)])
    ys=np.concatenate([np.full(nx+1,ylo),yb[1:],np.full(nx,1.0),yb[::-1][1:-1]])
    # exact maximal spacing (all coordinates are floats; distances bounded with outward rounding)
    dxs=up(np.abs(np.roll(xs,-1)-xs)); dys=up(np.abs(np.roll(ys,-1)-ys))
    hh=float(up(np.max(up(np.sqrt(up(dxs*dxs+dys*dys))))))
    return xs,ys,hh
# ------------------------------------------------------------------ time stepping
T0HI=hi(T0)                        # float >= t0: the last step ends at or beyond t0
def ylo_of(t2):                    # float <= sqrt(y0^2+2(t0-t)) for every t in [0, min(t2,t0)]
    if t2>=T0HI: return lo(Y0)
    return lo((Y0*Y0+2*(T0-arb(t2))).sqrt())
t=0.0; dt=2e-4; step=0; worst=1e9; npts=0; T=time.time(); h=None
while t<T0HI:
    t2=min(T0HI,t+dt)
    ylo=ylo_of(t2)
    Dz,Dt,E=bounds(t,t2,ylo)
    if h is None: h=min(0.02,1.0/Dz)
    xs,ys,hh=mesh(ylo,h)
    f,err=f_eval(xs,ys,t); m=float(dn(absmin(f)-err).min())
    if m<=Dz*hh:                      # mesh too coarse for a valid winding argument
        h=0.4*max(m,1e-6)/Dz; continue
    marg=m-Dz*hh/2-Dt*(t2-t)-E
    if marg<=0:
        dt=dt/2; assert dt>1e-10, 'cannot advance'; continue
    # rigorous winding number: (a) Re(f_{j+1} conj f_j) > 0 for every pair (interval arithmetic) => each true
    # increment of arg f lies in (-pi/2, pi/2); (b) the increments computed from the rectangle centres differ from
    # the true ones by <= asin(r_j/|f_j|)+asin(r_{j+1}/|f_{j+1}|) <= (pi/2)(r_j/m_j + r_{j+1}/m_{j+1}); the winding
    # number is an integer within (computed sum + total error) of 0, hence 0.
    g_=tuple(np.roll(v,-1) for v in f)
    rr=imul(g_[0],g_[1],f[0],f[1]); ii=imul(g_[2],g_[3],f[2],f[3])
    re_lo=dn(rr[0]+ii[0]); assert np.all(re_lo>0), 'consecutive values not within a quarter turn'
    mid=(f[0]+f[1])/2+1j*(f[2]+f[3])/2; inc=np.angle(np.roll(mid,-1)/mid)
    rad=up(np.sqrt(up(up((f[1]-f[0])/2)**2+up((f[3]-f[2])/2)**2))); mj=dn(absmin(f))
    werr=float(up(np.sum(up(2*(math.pi/2)*up(rad/mj)))*1.01+1e-12))
    wind=abs(float(np.sum(inc)))*(1+1e-12)+1e-6*len(inc)   # 1e-6 per term covers any libm error of np.angle
    assert wind+werr<math.pi, 'winding number not certified to be 0'
    worst=min(worst,marg); npts+=len(xs)
    print('step %3d t=[%.6f,%.6f] ylo=%.4f pts=%5d h=%.2e min|f|>=%.4f Dz<=%.1f Dt<=%.1f E<=%.1e margin=%.4f wind=0 [%.0fs]'%(step,t,t2,ylo,len(xs),hh,m,Dz,Dt,E,marg,time.time()-T),flush=True)
    t=t2; step+=1
    h=min(0.02,0.4*m/Dz)
    dt=max(1e-6,0.8*(m-Dz*h/2-E)/Dt)
assert t>=T0HI
print('CERTIFIED (iii) [Arb]: barrier X=%s, t in [0,%s] (last step ends at %.17g >= t0), y0=%s: no zeros of H_t; %d steps, %d evaluations, min margin %.4f, %.0fs'%(Xs,T0s,t,Y0s,step,npts,worst,time.time()-T))