Re: Truncation Precision in PARI

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

* Karim Belabas [2018-07-09 19:06]:
> * Karim Belabas [2018-07-09 18:32]:
> > * Peter Pein [2018-07-09 14:17]:
> > > that is not very accurate :(
> > > 
> > > I get:
> > > 
> > > precision(-1 - Euler + 3/2*log(2*Pi) + 6*zetahurwitz(-1, 1, der = 1), 100)
> > > %1 =
> > > 0.1870730725097797894509591576777666319578148029622159376465535484192711630046534855901322306210633101
[...]
> I misread your message, sorry : I had somehow entirely forgotten the
> original thread / question and the answer I made. Indeed,
> 
>   sumnum(k = 1, (1 - 3*k - 6*k^2) / (2*k+2) + 3*k^2 * log(1 + 1/k))
> 
> is very inaccurate. I only intended to show how intnum's expression
> could be transformed algbraically, it wasn't meant to be numerically
> stable or efficient and I didn't bother to check the numerical answer.
> As you pointed out, it's actually completely wrong. I'll investigate... 

The reason is "clear": we are summing

  f(k) = (1 - 3*k - 6*k^2) / (2*k+2) + 3*k^2 * log(1 + 1/k);

which is very inaccurate for large k due to catastrophic cancellation:
we are summing two terms of the order of O(k) for a result of the order
of O(k^(-2)). The sumnum() function calls intnum(), which itself uses the
double-exponential method and must evaluate f at very large values
indeed:

  ? \p100
  ? sumnum(k = 1, print(k); f(k))

will print values of the order of k ~ 5e136, for which no correct digit
is left after cancellation:
 
  ? f(1e136)
  %3 = 0.E22

(Other summation routines use much smaller values of k: ~33203 for
Monien summation and ~130 for Lagrange summation at \p100 accuracy.)

A brute-force fix is thus to replace the function evaluation by an
accurate asymptotic expansion for large k:
  ? F = truncate(f(1 / x + O(x^20)));
  ? sumnum(k = 1, if (k < 1e6,  f(k), subst(F,x,1/k)))
  %6 = 0.1870730725097797894509591576777666319578148029622159376465535484192711630046534855901322306210633101

A better fix would be to rewrite the summand so as to avoid cancellation...

Cheers,

    K.B.
--
Karim Belabas, IMB (UMR 5251)  Tel: (+33) (0)5 40 00 26 17
Universite de Bordeaux         Fax: (+33) (0)5 40 00 21 23
351, cours de la Liberation    http://www.math.u-bordeaux.fr/~kbelabas/
F-33405 Talence (France)       http://pari.math.u-bordeaux.fr/  [PARI/GP]
`