KEYWORDS: 
Binet function, Gamma Function, Continued Fraction, qd algorithm, 
Akiyama-Tanigawa algorithm, Bernoulli numbers, Jacques Binet, 
Thomas Jan Stieltjes, Heinz Rutishauser, Peter Henrici, 
OEIS A005146, OEIS A005147.

On Stieltjes' Continued Fraction for the Gamma Function.

[Math. of Comp., 34 (150), April 1980, 547-551]

This is the title of a paper of Bruce W. Char and is the
subject of this note.

J(z) = ln Gamma(z) + z - (z-1/2)ln(z) - ln sqrt(2 Pi)

J(z) is the Binet function. Stieltjes (1894) constructed
a continued fraction for J(z) and computed the first five
coefficients. He noted: "Le calcul est trés pénible ...
la loi de ces nombres paraît extrêmement compliqué".

Char used a recurrence relation based on an algorithm
of H.S. Wall to compute the first 40 coefficients
to 40 significant digits.

Peter Henrici used the qd-algorithm [Mathematical Methods
for Digital Computers, Vol.2, eds. Ralston and Wilf,
Wiley 1967] to sum asymptotic developments and gave a
Fortran program. This program is quite complicated and
I did not succeed to reconstruct it because of
out-of-range errors.

However, the following simple implementation of the
qd-algorithm with Maple is based on the same idea.

----------------------------------------------------
 qd := proc(S) local LEN,M,N,K,A,B,C;

 LEN := nops(S);
 M := array(1..LEN,1..LEN):

 for N to LEN   do M[N,1] := 0 od;
 for N to LEN-1 do M[N,2] := S[N+1]/S[N] od;

 for K from 3 to LEN do
   for N to LEN-K+1 do
     A := M[N+1,K-2]; B := M[N+1,K-1]; C := M[N,K-1];
     M[N,K] := `if`(type(K,odd), A+B-C, A*B/C);
   od;
 od;

 for N to LEN do M[N,1] := S[N] od;

 #print(M);
 [seq(M[1,K],K=1..LEN)] end:
----------------------------------------------------

Char used in 1979 the MACSYMA program and wrote:

"The computation of the coefficients a[0] through a[40],
as well as their 40-digit approximations, took approximately
twenty minutes of central processing unit time."

The same task can be performed with Maple and the above
routine as follows:

s := [seq((-1)^p*bernoulli(2*p+2)/((2*p+1)*(2*p+2)),p=0..40)];
t := time(); cf := qd(s); evalf(%,40); time()-t;

It took Maple V R5 only 8 seconds to accomplish this task.

==============================================================

It is interesting to compare this result with the 
following remark in:

"A Continued Fraction Algorithm for the Computation of 
Higher Transcendental Functions in the Complex Plane." 
I. Gargantini; P. Henrici, Mathematics of Computation,
Vol. 21, No. 97 (Jan., 1967), 18-29. 

5.1. Rational Arithmetic. If the coefficients c_n are 
rational, all entries of the QD-scheme are rational, 
and it is theoretically possible to construct the 
scheme in rational arithmetic by treating every number
as an ordered pair of integers [..]. This technique 
was found impractical, even in a simple case [..], 
because the number of digits in both numerator and 
denominator increases very rapidly with increasing k,
although every fraction was reduced to lowest terms 
by means of the Euclidian algorithm. Using a maximum
word length of 50 decimal digits on the IBM 1620 it 
was possible to compute the QD-scheme for K_0(z) only
up to k = 4. In another experiment, the QD-scheme 
associated with the function 

J(z) = log Gamma(z) — (z — 1/2) log z + z — (1/2) log 2 Pi 

(where the c_n are related to the Bernoulli numbers),
the calculations could only be carried out up to k = 3.

==============================================================

  C. Brezinski wrote in "Extrapolation algorithms and Padé
  approximations: a historical survey." (1994)
  "Related to the recursive computation of Padé approximants,
  is also the QD algorithm. It is due to Heinz Rutishauser
  but it was already contained is some earlier work by
  Stieltjes dating 1889. This algorithm was developed
  Peter Henrici."

==============================================================

If we look at the matrix M for LEN=5 we see the structure
of the computation. The first column is the input, the first
row the output of the algorithm. Viewed this way, we have a
general sequence to sequence transformation, valid for all
sequences s with s[n] <> 0 for all n.


     [  1         1         53         195       22999 ]
     [ ---       --         ---        ---       ----- ]
     [ 12        30         210        371       22737 ]
     [                                                 ]
     [  1         2         13        1841             ]
     [ ---        -         --        ----             ]
     [ 360        7         28        1716             ]
     [                                                 ]
     [  1         3         263                        ]
     [ ----       -         ---                        ]
     [ 1260       4         396                        ]
     [                                                 ]
     [  1        140                                   ]
     [ ----      ---                                   ]
     [ 1680       99                                   ]
     [                                                 ]
     [  1                                              ]
     [ ----                                            ]
     [ 1188                                            ]

In fact this matrix shows all what is to be done to arrive
at Stieltjes' five coefficients. Would he have known this
algorithm he certainly would not have written "Le calcul est
trés pénible ... " ;-).

============================================================

At the newsgroup sci.math.symbolic Richard J. Fateman wrote:

"Since the program does not require any special features of
computer algebra other than calculation of bernoulli numbers,
any system which has exact rational numbers, including any
Common Lisp implementation, should do just fine."

This remark inspired me to make the implementation independent
from the call to the Maple routine for the Bernoulli numbers.
So I added the Akiyama-Tanigawa algorithm for the calculation
of the Bernoulli numbers.

In this more self-contained form the procedure now is:

==============================================================

 StieltjesCF := proc(len) local s,m,n,k,a,b,c,AkiyamaTanigawa;

  # Computes Bernoulli(2*p+2)/((2*p+1)*(2*p+2))
  # using the Akiyama-Tanigawa algorithm

  AkiyamaTanigawa := proc(l) local a,t,n,m,j,k;
  n := 2*l+1;
  t := array(1..n);
  a := array(1..l);
  t[1] := 1; k := 1;
  for m from 2 to n do
    t[m] := 1/m;
    for j from m-1 by -1 to 1 do
      t[j] := j*(t[j]-t[j+1]);
    od;
    if type(m,odd) then
       a[k] := (-1)^(k+1)*t[1]/((2*k-1)*(2*k));
       k := k+1;
    fi;
  od;
  a end:

  # Computes the Stieltjes continued fraction for the
  # Gamma function using Rutishauser's QD-algorithm.

  s := AkiyamaTanigawa(len);
  m := array(1..len,1..len):
  for n to len   do m[n,1] := 0 od;
  for n to len-1 do m[n,2] := s[n+1]/s[n] od;
  for k from 3 to len do
    for n to len-k+1 do
      a := m[n+1,k-2]; b := m[n+1,k-1]; c := m[n,k-1];
      m[n,k] := `if`(type(k,odd), a+b-c, a*b/c);
    od;
  od;
  m[1,1] := s[1];
  seq(m[1,k],k=1..len) end:

  # Example call
  StieltjesCF(6);

      1   1  53   195  22999  29944523
     --, --, ---, ---, -----, --------
     12  30  210  371  22737  19733142

=============================================================
See also at the "On-Line Encyclopedia of Integer Sequences":

http://www.research.att.com/~njas/sequences/?q=A005146
http://www.research.att.com/~njas/sequences/?q=A005147

A005146 := proc(n) numer(StieltjesCF(n)) end;
A005147 := proc(n) denom(StieltjesCF(n)) end;

==============================================================

This appendix shows how much more readable a Maple or Maxima 
program is compared to an implementation in a conventional
mainstream programming language like C#. If implemented in
Java the program turns almost unreadable because Java does
not support operator overloading. Moreover, it was 
interesting to observe that the compiled version of the 
following program was not faster than its Maple counterpart.

Appendix: An implementation in Csharp
=====================================
       
using System;
using Rational;
using RatMath;

namespace ContFrac
{
    public class ContinuedFraction
    {
        // Computes the sequence of Bernoulli(2*k) numbers (k >= 1)
        // using the Akiyama-Tanigawa algorithm

        static Rational[] bernoulli(int l)
        {
            int n = 2 * l + 1, k = 0; 
            Rational[] t = new Rational[n];
            Rational[] a = new Rational[l];

            t[0] = RatMath.of_int(1);

            for (int m = 1; m < n; m++)
            {
                t[m] = RatMath.rational(1, m + 1);
                for (int j = m; j > 0; j--)
                {
                    t[j - 1] = RatMath.of_int(j) * (t[j - 1] - t[j]);
                }

                // Get all Bernoulli numbers by deleting the 'if' clause.
                if ((m & 1) == 0) { a[k++] = t[0]; }
            }
            return a;
        }

        // Computes Rutishauser's Quotient-Difference (QD-) algorithm

        static Rational[] qd(Rational[] s)
        {
            int len = s.Length;
            Rational a, b, c;
            Rational zero = RatMath.of_int(0);
            Rational[,] m = new Rational[len, len];
            Rational[] r = new Rational[len];
            
            for (int n = 0; n < len; n++) { m[n, 0] = zero; }
            for (int n = 0; n < len - 1; n++) { m[n, 1] = s[n+1] / s[n]; }

            r[0] = s[0]; r[1] = m[0, 1];

            for (int k = 2; k < len; k++)
            {
                for (int n = 0; n < len - k ; n++)
                {
                    a = m[n + 1, k - 2]; b = m[n + 1, k - 1]; c = m[n, k - 1];
                    m[n, k] = (k & 1) == 0 ? a + b - c : a * b / c;
                }
                r[k] = m[0, k];
            }
            return r;
        }

        // Decorates the even Bernoulli numbers with weights (-1)^k/((2*k+1)*(2*k+2))

        static Rational[] weights(Rational[] s)
        {
            int sgn = 1;
            for (int k = 0; k < s.Length; k++)
            {
                s[k] = s[k] * RatMath.rational(sgn, (2 * k + 1) * (2 * k + 2));
                sgn = -sgn;
            }
            return s;
        }

        // Computes the Stieltjes continued fraction for the
        // Gamma function using Rutishauser's QD-algorithm.

        static Rational[] StieltjesCF(int len)
        {
            return qd(weights(bernoulli(len)));
        }

        static void Main(string[] args)
        {
            int len = 7;
            Rational[] cf = StieltjesCF(len);
            for (int k = 0; k < len; k++) 
               Console.WriteLine(k + ".." + cf[k]); 
            Console.ReadLine();
       }
    }
}       

previous Approximations for the Factorial Function.