from mpmath import mp, cos, sin, exp, log, pi, e, findroot
mp.dps = 50


#tc = 49.7738324776
tc = 100.0

def hardy_Z(t_val, method="siegelz"):
    """
    Hardy Z(t).

    `method="siegelz"` uses mpmath's stable implementation.
    `method="integral"` uses the eta-integral representation for validation,
    but it is much less stable for larger t.
    """
    t_val = mp.mpf(t_val)
    if method == "integral":
        s = mp.mpf('0.5') + 1j * t_val
        zeta_s = zeta_via_integral(s)
        theta = mp.siegeltheta(t_val)
        return mp.re(mp.e**(1j * theta) * zeta_s)
    return mp.siegelz(t_val)

def Z_func(t_val):
    """Stable Hardy Z(t) used for root finding."""
    return hardy_Z(t_val)
    
t_current = tc
for _ in range(10):
    Z_val = Z_func(t_current)
    # Numerical derivative
    dt = mp.mpf('1e-10')
    Z_prime = (Z_func(t_current + dt) - Z_func(t_current - dt)) / (2 * dt)
    
    if abs(Z_prime) < 1e-100:
        break
    
    t_next = t_current - Z_val / Z_prime
    if abs(t_next - t_current) < 1e-50:
        print(t_current, Z_val)
        break
    t_current = t_next