Gamma and Bernoulli

[Dirk Laurie <dpl@sun.ac.za>]

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).

I tried the following experiment:
? allocatemem(400000000)
? \p1000
   realprecision = 1001 significant digits (1000 digits displayed)
? #
   timer = 1 (on)
? gamma(0.785);
time = 10,999 ms.
? gamma(0.786);
time = 214 ms.
? bernreal(2000)
time = 0 ms.
? bernreal(3000);
time = 3,214 ms.
? bernreal(2998);
time = 0 ms.
? bernfrac(2000);
time = 1,678 ms.
? bernfrac(2000);
time = 1,688 ms.

I deduce the following:
  1. Yes, Pari-GP uses the asymptotic expansion.
  2. It remembers the real-valued Bernoulli numbers, so the second
     value of gamma is much faster than the first.
  3. It requires between 2000 and 3000 Bernoulli numbers at precision
     1000.  
  4. It uses the recursive algorithm to calculate those numbers,
     since it can deliver bernreal(2998) instantly once bernreal(3000)
     is known.
  5. It does not remember rational Bernoulli numbers.

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!

Dirk