Re: Gamma and Bernoulli

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

On Thu, 6 Nov 2003, Dirk Laurie wrote:
> I have an application in which I often need very accurate (about
> 10000 digits) of gamma.  The obvious algorithm seems to be the
> aymptotic expansion of y=log(gamma(x)) in terms of Bernoulli numbers
> coupled with moving the argument to a fairly large value so the
> expansion converges fast, recursing back, and evaluating exp(y).

Indeed that's the algorithm used.

>   4. It uses the recursive algorithm to calculate those numbers,
>      since it can deliver bernreal(2998) instantly once bernreal(3000)
>      is known.

Not anymore. In the the CVS version, a few routines compute individual
Bernoulli numbers of high index k using the Euler product for zeta(2*k).

>   5. It does not remember rational Bernoulli numbers.

Which is a bit silly. A more efficient algorithm would be to cache floating
point Bernoulli numbers B_2i, as is done currently, but also cache it in
exact rational form as soon as precision becomes large enough to recognize
B_2i ( since its denominator is known and small ). This way

1) increasing cache precision would be rather painless (one division per
number, and the denominator is small).

2) computations in large accuracy could use these rational numbers of small
height for most of the computation, and switch to floats only when the
height of B_2i becomes larger than the relative accuracy.


Food for thought: to compute a huge number of Bernoullis, a Newton inversion
(using FFT multiplication) of the truncated power series for (exp(x) - 1) / x
could be more efficient than the standard iteration.

> I would very much like to save time by saving Bernoulli numbers so
> I can just read them in from a file next time.  The computation of
> that first gamma value at 10000 digits is a bit time-consuming!

This is a nice idea, but currently not possible. Something like

  bernset( [1/6, -1/30, ... ] );

[ or floating point approximations ] would be nice indeed.

I'll include this in the TODO list.

Ceers,

    Karim.

P.S: the gamma() function should be quite a bit faster in the CVS version
( and require less Bernoulli ).
-- 
Karim Belabas                     Tel: (+33) (0)1 69 15 57 48
Dép. de Mathématiques, Bât. 425   Fax: (+33) (0)1 69 15 60 19
Université Paris-Sud              http://www.math.u-psud.fr/~belabas/
F-91405 Orsay (France)            http://pari.math.u-bordeaux.fr/  [PARI/GP]