Re: Question on finding a Riemann Zeta function zero for high values of s

[Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>]

On Thu, Oct 23, 2025 at 09:15:45AM +0200, Cohen Henri wrote:
> 
> Sorry, you are right, my patch breaks things. The unpatched code does give the correct
> result. I will look into it again.

I suggest a different patch below where instead of
exp(p1)*p2 we compute exp(p1+log(p2))
so cancellation between p1 and log(p2) can occur.

See the branch henri-lfunlarge_overflow

Cheers,
Bill.

diff --git a/src/basemath/lfunlarge.c b/src/basemath/lfunlarge.c
index 686738a240..91f9524553 100644
--- a/src/basemath/lfunlarge.c
+++ b/src/basemath/lfunlarge.c
@@ -301,11 +301,12 @@ integrand_h0(GEN sel, GEN s, GEN VCALL, GEN x, long prec)
 {
   pari_sp av = avma;
   long md = get_modulus(VCALL);
-  GEN r0 = m_r0(sel), aleps = m_aleps(sel), zn, p1;
+  GEN r0 = m_r0(sel), aleps = m_aleps(sel), zn, p1, p2;
   GEN pmd = divru(mppi(prec), md), ix = ginv(x);
   zn = gadd(r0, gdivgs(gmul(aleps, gsub(x, ix)), 2));
-  p1 = gmul(expIxy(pmd, gsqr(zn), prec),
-            gmul(gpow(zn, gneg(s), prec), gmul(aleps, gadd(x, ix))));
+  p1 = gmul(pmd, mulcxI(gsqr(zn)));
+  p2 = gadd(gmul(glog(zn,prec), gneg(s)), gadd(glog(aleps,prec), glog(gadd(x, ix),prec)));
+  p1 = gexp(gadd(p1,p2), prec);
   if (md == 1)
     p1 = gdiv(mkvec(mulcxI(p1)), gmul2n(gsin(gmul(pmd, zn), prec), 2));
   else