import mpmath as mm
mm.mp.dps=40
NS=300000; sv=bytearray([1])*(NS+1); sv[0]=sv[1]=0
for i in range(2,int(NS**.5)+1):
    if sv[i]: sv[i*i::i]=bytearray(len(sv[i*i::i]))
P=[p for p in range(NS+1) if sv[p]]
rows=[l.split() for l in open('../sieve/cp_1e9_1e10.txt')]
worst=-1e9
for r in rows[::4]:
    xi=int(r[0]); th=mm.mpf(r[1]); M=mm.mpf(r[2]); Q2=mm.mpf(r[3])
    t=mm.log(xi); eps=mm.log1p(mm.mpf(1)/xi)/t
    def xik(k):
        f=lambda u: mm.log1p(1/(mm.e**u*(mm.e**(k*u)-1)/(mm.e**u-1)))-eps*u
        return mm.e**mm.findroot(f,(mm.log(1.5),(t+mm.log(k))/k+1),solver='anderson')
    xs={};k=2
    while True:
        if mm.log1p(1/(2*(2**k-1)))/mm.log(2) < eps: break
        z=xik(k)
        xs[k]=z;k+=1
    a={p:max(k for k in xs if p<=xs[k]) for p in P if p<=xs[2]}
    logN=th+sum(sum(mm.log(p) for p in P if p<=xs[k]) for k in xs)
    Lam=(Q2-sum(mm.log(1-mm.mpf(p)**-2) for p in P if p<=xs[2]))+sum(mm.log(1-mm.mpf(p)**(-a[p]-1)) for p in a)
    ls=-M+Lam; g=mm.euler
    c=(mm.e**(ls-g)-mm.log(logN))*mm.sqrt(logN)
    a2=a[2]; logM=logN-a2*mm.log(2); lsM=ls-mm.log(2*(1-mm.mpf(2)**(-a2-1)))
    codd=(2*mm.e**(lsM-g)-mm.log(logM))*mm.sqrt(logM)
    off=(codd-c)*mm.sqrt(xi)
    bound=t+mm.log(t)+0.37
    worst=max(worst,float(off-bound))
    print(f"{xi:.4e} a={a2} offset*sqrt(xi)={float(off):.4f} bound_0.37={float(bound):.4f} ratio_to_t={float(off/t):.4f} {'VIOLATED' if off>bound else ''}")
print("max(offset - bound)*sqrt(xi):",worst)