Stieltjes Continued Fraction for the Gamma Function
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(); } } } 
Source: http://www.luschny.de/math/factorial/approx/continuedfraction.html
0 Response to "Stieltjes Continued Fraction for the Gamma Function"
Post a Comment