"""Heat flow = Gaussian average of shifted Dirichlet sums:  e^{(t/4) L^2} = E[e^{cZL}], c = sqrt(t/2), Z ~ N(0,1).
So  A = sum_{n<=N} b_n n^{-s*} = E_Z[ X_N(s*-cZ) ],  X_N partial zeta sum (on the phase torus).
For sigma' = Re(s*)-cZ >= sigma_c > 1, X_N = (Euler product) + tail, and on the whole torus
   Re zeta_theta(sigma') >= zeta(2 sigma')/zeta(sigma') * cos Phi(sigma'),  Phi(sigma) = sum_p arcsin(p^-sigma) < pi/2.
Lower bound for Re A (all phases at once) and upper bound for |gamma||B|; compare with the triangle / mollified bounds."""
import numpy as np, sys
from scipy.special import zeta
from scipy.stats import norm
def primes(n):
    s=np.ones(n+1,bool); s[:2]=False
    for p in range(2,int(n**0.5)+1):
        if s[p]: s[p*p::p]=False
    return np.nonzero(s)[0].astype(float)
PR=primes(2_000_000)
def Phi(sig):     # sum_p arcsin(p^-sig) (prime tail beyond 2e6 bounded by integral)
    v=np.sum(np.arcsin(PR**(-sig)))
    P0=PR[-1]; v+= P0**(1-sig)/((sig-1)*np.log(P0))*1.1
    return v
def bounds(x,y,t,sig_c=1.25,N=None):
    s_star=(1+y)/2+(t/4)*np.log(x/(4*np.pi)); c=np.sqrt(t/2)
    if N is None: N=int(np.sqrt(x/(4*np.pi)))
    z=np.linspace(-8,8,16001); w=norm.pdf(z)*(z[1]-z[0])
    sp=s_star-c*z                                   # sigma' for the A part
    good=sp>=sig_c
    lowA=np.zeros_like(z); badA=np.zeros_like(z)
    for i in np.nonzero(good)[0]:
        sg=sp[i]; ph=Phi(sg)
        if ph>=np.pi/2: good[i]=False; continue
        tail=N**(1-sg)/(sg-1)
        lowA[i]=zeta(2*sg)/zeta(sg)*np.cos(ph)-tail
    n=np.arange(1,N+1,dtype=float); ln=np.log(n)
    for i in np.nonzero(~good)[0]:
        badA[i]=np.sum(np.exp(-sp[i]*ln))               # |X_N| <= sum n^-sigma'
    ReA_lo=np.sum(w*lowA)-np.sum(w*badA)
    # B part: |gamma| |sum n^y b_n n^{-s*-...}| <= |gamma| E_Z[ sum n^{-(sigma*-y-cZ)} ] ; use zeta(sigma'') where sigma''>1+eps
    g=np.exp(0.02*y)*(x/(4*np.pi))**(-y/2)
    spB=s_star-y-c*z; upB=np.zeros_like(z)
    for i,sg in enumerate(spB):
        upB[i]= (zeta(sg)+0.0) if sg>1.05 else np.sum(np.exp(-sg*ln))
    B_up=g*np.sum(w*upB)
    return s_star,c,ReA_lo,np.sum(w[~good]),B_up,ReA_lo-B_up
if __name__=='__main__':
    x=float(sys.argv[1]); y=float(sys.argv[2])
    for t in [float(v) for v in sys.argv[3:]]:
        for sc in [1.15,1.2,1.25,1.3,1.4]:
            s_star,c,RA,wb,BU,LB=bounds(x,y,t,sc)
            print('x=%.1e y=%.2f t=%.3f sigma*=%.3f c=%.3f sigma_c=%.2f: ReA>=%.4f (bad weight %.3f)  |gB|<=%.4f   LB=%+.4f'%(x,y,t,s_star,c,sc,RA,wb,BU,LB),flush=True)