Re: Euler product

[Karim Belabas <Karim.Belabas@math.u-psud.fr>]

* old [2004-05-03 12:08]:
> Here is the constants I want to compute:
> 
> C(h)=prodeuler(p=2,infinity, (1+(1+1/p)^h/(p-1))/(1+(1+1/p)^(h-1)/(p-1)))
> 
> and with h quite large (up to 50).
> 
> I have a prodeulerrat by Henri Cohen, but the results are incoherent

[ which may be downloaded from my GP scripts page as

    http://www.math.u-psud.fr/~belabas/pari/scripts/cohen.gp
]

> from h=19-21 onwards. I tried to compute C(h)/C(h-1) directly but
> failed also. Log is below.
[...]
> {frel(h)=prodeulerrat((1+(p+1)^h/p^h/(p-1))/(1+(p+1)^(h-1)/p^(h-1)/(p-1)),2)}
> 
> {g(n)=local(res,resrel,h); res=1.519817754635066571658191948;
>            for(h=2,n,resrel=frel(h);res*=resrel;print(h," -> 
> ",resrel,"  --> ", res));}
> 
> ? g(50)
> 2 -> 1.628355734896019122  --> 2.474803956756801499
[...]
> 21 -> 2.797113371058406802  --> 35081822.12856861096
> 22 -> 2.423148937773094341  --> 85008480.02598566512
> 23 -> 1.541937578649181143  --> 131077769.8559156157
> 24 -> 0.003825393288703888199  --> 501424.0211050924189
> 25 -> 5.443317760119889255 E32  --> 2.729410279432079743 E38

Looks like prodeulerrat uses heuristic bounds to evaluate the required
accuracy. Catastrophic cancellation occured and was hidden away by the
subsequent precision(res, pr), which added meaninless trailing zeroes to
the result.

Try the attached variant, it should be slightly slower, but much safer
[ also includes sumeulerrat() ].

Cheers,

    Karim.
-- 
Karim Belabas                     Tel: (+33) (0)1 69 15 57 48
Dep. de Mathematiques, Bat. 425   Fax: (+33) (0)1 69 15 60 19
Universite Paris-Sud              http://www.math.u-psud.fr/~belabas/ 
F-91405 Orsay (France)            http://pari.math.u-bordeaux.fr/  [PARI/GP]

logzetaa(N,A) = log(zeta(N) * prodeuler(p=2,A,1-1/p^N))
eval_sal(sal, A, pr) =
{ local(v, m); \\ for efficiency, realprecision had better be low at this point

  v = abs( Vec(sal * 1.) ); m = 0;
  for (i = 1, length(v), if (v[i] && v[i] < m, m = v[i]));
  if (!m, m = 1);
  default(realprecision,
          pr + max(0, floor(log(vecmax(v) / m)/log(10))));
  sal *= 1.; 
  sum(N = 2, lim, logzetaa(N,A) * sumdiv(N,k,moebius(k)*polcoeff(sal,N/k)/k));
} 
restore_pr(res, pr) =
{ default(realprecision, pr);
  if (precision(res) > pr, precision(res,pr), res);
} 

/* F rational function. Computes the sum over primes from p=pin to infinity
   of F(p). We must have F(x)=O(1/x^2) as x tends to infty. */
sumeulerrat(F, pin, A = 30) =
{ local(vx, pr, mro, lim, sal, v, res, m);

  vx = variable(F); pr = precision(1.);
  if (poldegree(F) > -2, error("infinite sum in sumeulerrat"));
  default(realprecision, 9);
  mro = max(1, vecmax(abs(polroots(denominator(F)))));
  A = max(pin, ceil(max(A, 3*mro)));
  lim = 2 + ceil(pr*log(10) / log(A/mro));
  sal = subst(F,vx,1/x) + O(x^lim); 
  res = eval_sal(sal, A, pr);
  forprime (p = pin, A, res += subst(F,vx,p));
  restore_pr(res, pr);
} 
  
/* F rational function. Computes the product over primes from p=pin to infinity
   of F(p). We must have F(x)=1+O(1/x^2) as x tends to infty. */

prodeulerrat(F, pin, A = 30) =
{ local(vx, pr, mro, lim, sal, v, res);

  vx = variable(F); pr = precision(1.);
  if (poldegree(F-1) > -2, error("infinite product in prodrat"));
  default(realprecision, 9);
  mro = max(1, vecmax(abs(polroots(denominator(F)))));
  mro = max(mro, vecmax(abs(polroots(numerator(F)))));
  A = max(pin, ceil(max(A, 3*mro)));
  lim = 2 + ceil(pr*log(10) / log(A/mro));
  sal = log(subst(F,vx,1/x) + O(x^lim)); 
  res = exp( eval_sal(sal, A, pr) );
  forprime (p = pin, A, res *= subst(F,vx,p));
  restore_pr(res, pr);
}