AN ANALYSIS OF THE LANCZOS GAMMA

APPROXIMATION
by

GLENDON RALPH PUGH

B.Sc. (Mathematics) University of New Brunswick, 1992

M.Sc. (Mathematics) University of British Columbia, 1999

A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF

THE REQUIREMENTS FOR THE DEGREE OF

DOCTOR OF PHILOSOPHY

in

THE FACULTY OF GRADUATE STUDIES

Department of Mathematics

We accept this thesis as conforming

to the required standard

..................................

..................................

..................................

..................................

..................................

THE UNIVERSITY OF BRITISH COLUMBIA

November 2004

c Glendon Ralph Pugh, 2004

In presenting this thesis in partial fulfillment of the requirements for an ad-
vanced degree at the University of British Columbia, I agree that the Library
shall make it freely available for reference and study. I further agree that
permission for extensive copying of this thesis for scholarly purposes may be
granted by the head of my department or by his or her representatives. It is
understood that copying or publication of this thesis for financial gain shall
not be allowed without my written permission.

(Signature)

Department of Mathematics
The University of British Columbia
Vancouver, Canada

Date
Abstract
This thesis is an analysis of C. Lanczos’ approximation of the classical
gamma function Γ(z +1) as given in his 1964 paper A Precision Approximation
of the Gamma Function [14]. The purposes of this study are:

(i) to explain the details of Lanczos’ paper, including proofs of all claims
made by the author;
(ii) to address the question of how best to implement the approximation
method in practice; and
(iii) to generalize the methods used in the derivation of the approximation.

At present there are a number of algorithms for approximating the gamma

function. The oldest and most well-known is Stirling’s asymptotic series which
is still widely used today. Another more recent method is that of Spouge [27],
which is similar in form though different in origin than Lanczos’ formula. All
three of these approximation methods take the form
√
Γ(z + 1) = 2π(z + w)z+1/2 e−z−w [sw,n (z) + ǫw,n (z)] (1)

where sw,n (z) denotes a series of n + 1 terms and ǫw,n (z) a relative error to
be estimated. The real variable w is a free parameter which can be adjusted
to control the accuracy of the approximation. Lanczos’ method stands apart
from the other two in that, with w ≥ 0 fixed, as n → ∞ the series sw,n (z)
converges while ǫw,n (z) → 0 uniformly on Re(z) > −w. Stirling’s and Spouge’s
methods do not share this property.
What is new here is a simple empirical method for bounding the relative
error |ǫw,n (z)| in the right half plane based on the behaviour of this function
as |z| → ∞. This method is used to produce pairs (n, w) which give formu-
las (1) which, in the case of a uniformly bounded error, are more efficient than
Stirling’s and Spouge’s methods at comparable accuracy. In the n = 0 case,
a variation of Stirling’s formula is given which has an empirically determined
uniform error bound of 0.006 on Re(z) ≥ 0. Another result is a proof of the
limit formula

z 1 −12 /r z −22 /r z(z − 1)
Γ(z + 1) = 2 lim r −e +e +···
r→∞ 2 z+1 (z + 1)(z + 2)

as stated without proof by Lanczos at the end of his paper.

ii
Table of Contents
Abstract ii

Table of Contents iii

List of Tables vi

List of Figures vii

Acknowledgments viii

Chapter 1. Introduction 1
1.1 Lanczos and His Formula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

Chapter 2. A Primer on the Gamma Function 20

2.1 First Encounters with Gamma . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.2 A Brief History of Gamma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.3 Definition of Γ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.4 Standard Results and Identities . . . . . . . . . . . . . . . . . . . . . . . . 27
2.5 Stirling’s Series and Formula . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.5.1 The Formula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.5.2 Some Remarks about Error Bounds . . . . . . . . . . . . . . 34
2.6 Spouge’s Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.7 Additional Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

Chapter 3. The Lanczos Formula 40

3.1 The Lanczos Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.1.1 Preliminary Transformations . . . . . . . . . . . . . . . . . . . . . 42
3.1.2 The Implicit Function v(θ) . . . . . . . . . . . . . . . . . . . . . . 44
3.1.3 Series Expansion of fE,r . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.2 Derivation in the Chebyshev Setting . . . . . . . . . . . . . . . . . . . . 47
3.3 Derivation in the Hilbert Space Setting . . . . . . . . . . . . . . . . . 49
3.4 Additional Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

Chapter 4. The Functions v and fE,r 52

4.1 Closed Form Formulas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.1.1 Closed Form for v(x) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

iii
Table of Contents

4.1.2 Closed Form for fr (x) and fE,r (x) . . . . . . . . . . . . . . . . 54

4.2 Smoothness of v r (x), fr (x) and fE,r (x) . . . . . . . . . . . . . . . . . 56
4.3 Smoothness of v r (θ), fr (θ) and fE,r (θ) . . . . . . . . . . . . . . . . . . 59
4.3.1 Generalized fr (θ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.3.2 Derivatives of G(β, j, k, l; θ) . . . . . . . . . . . . . . . . . . . . . . 61
4.3.3 Application to fr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.4 Fourier Series and Growth Order of Coefficients . . . . . . . . . 64

Chapter 5. Analytic Continuation of Main Formula 66

5.1 Domain of Convergence of Sr (z) . . . . . . . . . . . . . . . . . . . . . . . . 66
5.2 Analytic Continuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.3 Additional Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

Chapter 6. Computing the Coefficients 72

6.1 The Lanczos Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
6.2 A Horner Type Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
6.3 Another Recursion Formula . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.4 Discrete Fourier Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.4.1 Finite Fourier Series . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.4.2 Finite Chebyshev Series . . . . . . . . . . . . . . . . . . . . . . . . . 78
6.5 Some Sample Coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6.6 Partial Fraction Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . 80
6.7 Matrix Representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

Chapter 7. Lanczos’ Limit Formula 84

7.1 Some Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
7.2 Proof of the Limit Formula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
7.3 Additional Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

Chapter 8. Implementation & Error Discussion 98

8.1 A Word or Two about Errors . . . . . . . . . . . . . . . . . . . . . . . . . . 99
8.1.1 Absolute and Relative Error . . . . . . . . . . . . . . . . . . . . . 99
8.1.2 Lanczos’ Relative Error . . . . . . . . . . . . . . . . . . . . . . . . . 100
8.1.3 Restriction to Re(z) ≥ 0 . . . . . . . . . . . . . . . . . . . . . . . . 101
8.2 Some Notation and a Theorem . . . . . . . . . . . . . . . . . . . . . . . . 102
8.3 Lanczos’ Error Estimates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
8.4 Luke’s Error Estimates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
8.5 Numerical Investigations & the Zeros of ǫ∞ r,n . . . . . . . . . . . 109
8.5.1 The Motivating Idea . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
8.5.2 Zeros of ǫ∞
r,n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

iv
Table of Contents

8.5.3 Improved Lanczos Formula . . . . . . . . . . . . . . . . . . . . . 111

8.5.4 Optimal Formulas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
8.6 Additional Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

Chapter 9. Comparison of the Methods 120

9.1 Stirling’s Series . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
9.1.1 Γ(20 + 17i) with |ǫρ | < 10−32 . . . . . . . . . . . . . . . . . . . 120
9.1.2 Γ(1 + z) with |ǫρ | < 10−32 Uniformly . . . . . . . . . . . . 122
9.2 Spouge’s Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
9.2.1 Γ(20 + 17i) with |ǫρ | < 10−32 . . . . . . . . . . . . . . . . . . . 123
9.2.2 Γ(1 + z) with |ǫρ | < 10−32 Uniformly . . . . . . . . . . . . 124
9.3 Lanczos’ Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
9.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

Chapter 10. Consequences of Lanczos’ Paper 128

10.1 The Lanczos Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
10.2 A Combinatorial Identity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
10.3 Variations on Stirling’s Formula . . . . . . . . . . . . . . . . . . . . . . . 134

Chapter 11. Conclusion and Future Work 139

11.1 Outstanding Questions Related to Lanczos’ Paper . . . . . 140
11.2 Outstanding Questions Related to Error Bounds . . . . . . . 141
11.3 Conjectures on a Deterministic Algorithm . . . . . . . . . . . . . 143

Bibliography 145

Appendix A. Preliminary Steps of Lanczos Derivation 148

Appendix B. A Primer on Chebyshev Polynomials 149

Appendix C. Empirical Error Data 152

v
List of Tables
1.1 Relative decrease of ak (r), r = 1, 4, 7 . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.1 Upper bound on |EN,n (6 + 13i)| in Stirling Series . . . . . . . . . . . . 33

2.2 Coefficients of Spouge’s Series, a = 12.5 . . . . . . . . . . . . . . . . . . . . . 38

6.1 Scaled summands of c6 (6) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

6.2 Coefficients as functions of r . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

8.1 Lanczos’ uniform error bounds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

8.2 Location of maximum of ηr,n (θ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
8.3 Errors at infinity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
8.4 Approximate Zeros of ǫ∞ r,6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
8.5 Sample coefficients, n = 10, r = 10.900511 . . . . . . . . . . . . . . . . . 116
8.6 Comparison of Mr(n),n and an+1 (r(n)) . . . . . . . . . . . . . . . . . . . . . . 119

9.1 Upper bound on |EN,n (19 + 17i)| in Stirling Series . . . . . . . . . 121

9.2 Upper bound on |EN,n (0)| in Stirling Series . . . . . . . . . . . . . . . . 122
9.3 Coefficients of Spouge’s series, a = 39 . . . . . . . . . . . . . . . . . . . . . . 124
9.4 Coefficients of Lanczos’ series, r = 22.618910 . . . . . . . . . . . . . . . 126

C.1 Error Data, n = 0, . . . , 30 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153

C.2 Error Data, n = 31, . . . , 60 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154

vi
List of Figures
1.1 Cornelius Lanczos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Relative decrease of ak (1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Relative decrease of ak (4) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4 Relative decrease of ak (7) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.5 Lanczos shelf, r = 20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.6 Lanczos shelf, r = 40 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.7 Lanczos shelf, r = 60 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.8 fr (θ), −π/2 ≤ θ ≤ π/2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.1 Γ(z + 1), −4 < z < 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

2.2 |Γ(x + iy)|, −3 < x < 4, −3 < y < 3 . . . . . . . . . . . . . . . . . . . . . . . 28

3.1 Path of v(1 − log v) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

4.1 v(x)/e, −1 ≤ x ≤ 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
√
4.2 2 fE,r (x)/er , −1 ≤ x ≤ 1, r = 0, 1, 2, 4 . . . . . . . . . . . . . . . . . . . . 55

5.1 Growth of ak (r), k −(2r+1) , and k −⌈2r⌉ with r = 1.25 . . . . . . . . . 71

7.1 exp [−r(θ − π/2)2 ] − gr (θ) for r = 1, 6, 10 . . . . . . . . . . . . . . . . . . . 86

7.2 log Q(k, 10), k = 0, . . . , 24 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
7.3 log Q(k, 20), k = 0, . . . , 24 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
7.4 log Q(k, 50), k = 0, . . . , 24 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
∞ 1/5
8.1 ǫr,n , n = 0, 1, 2, 3, −1/2 < r < 6 . . . . . . . . . . . . . . . . . . . . . . 112
8.2 log ǫrj ,6 (it/(1 − t)), 0 < t < 1, j = 0, . . . , 11 . . . . . . . . . . . . . . 114

8.3 log |ǫr,4 (it/(1 − t))|, 0 < t < 1, 0 < r < 5 . . . . . . . . . . . . . . . . . . 115
8.4 [ǫ∞
r,6 ]
1/7
, −0.2 < r < 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
1/7
8.5 [a7 (r)] , −0.2 < r < 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
8.6 [ǫ∞
r,6 ]
1/7
and [a7 (r)]1/7 , −0.2 < r < 7 . . . . . . . . . . . . . . . . . . . . . . . 118
z z
10.1 Path of a(0, 1−z ), a(−1, 1−z ), 0 ≤ z < 1 . . . . . . . . . . . . . . . . . . . 138

11.1 ǫr,n (z), −0.5 < r < 3, −0.5 < z < 3 . . . . . . . . . . . . . . . . . . . . . . . 142
11.2 Linear fit of (n, − log Mr(n),n ), n = 0, . . . , 60 . . . . . . . . . . . . . . . 144

vii
Acknowledgments
I would like to thank all those who have assisted, guided and sup-
ported me in my studies leading to this thesis. In particular, I wish
to thank my supervisor Dr. Bill Casselman for his patience, guidance
and encouragement, my committee members Dr. David Boyd and Dr.
Richard Froese, and Ms. Lee Yupitun for her assistance as I found my
way through the administrative maze of graduate school. I would also
like to thank Professor Wesley Doggett of North Carolina State Uni-
versity for his kind assistance in obtaining permission to reproduce the
photo of C.Lanczos in Chapter 1.
The online computing community has been a great help to me also
in making available on the web many useful tools which helped in my
research. In particular, MPJAVA, a multiple precision floating point
computation package in Java produced by The HARPOON Project
at the University of North Carolina was used for many of the high
precision numerical investigations. In addition, the online spreadsheet
MathSheet developed by Dr. Bill Casselman and Dr. David Austin was
very helpful in producing many of the PostScript plots in the thesis.
Finally, I would like to thank my family for their support and en-
couragement: my parents Anne and Charles, and especially my wife
Jacqueline for her encouragement and many years of understanding
during the course of my studies.

viii
Chapter 1
Introduction

The focus of this study is the little known Lanczos approximation [14]
for computing the gamma function Γ(z + 1) with complex argument z.
The aim is to address the following questions:

(i) What are the properties of the approximation? This involves a

detailed examination of Lanczos’ derivation of the method.

(ii) How are the parameters of the method best chosen to approximate
Γ(z + 1)?

(iii) How does Lanczos’ method compare against other methods?

(iv) Does Lanczos’ method generalize to other functions?

The subject of approximating the factorial function and its gener-

alization, the gamma function, is a very old one. Within a year of
Euler’s first consideration of the gamma function in 1729, Stirling pub-
lished the asymptotic formula for the factorial with integer argument
which bears his name. The generalization of this idea, Stirling’s series,
has received much attention over the years, especially since the advent
of computers in the twentieth century.
Stirling’s series remains today the state of the art and forms the
basis of many, if not most algorithms for computing the gamma func-
tion, and is the subject of many papers dealing with optimal computa-
tional strategies. There are, however, other methods for computing the
gamma function with complex argument [7][20][26][14][27]. One such
method is that of Lanczos, and another similar in form though different
in origin is that of Spouge.

1
Chapter 1. Introduction

The methods of Stirling and Lanczos share the common strategy

of terminating an infinite series and estimating the error which results.
Also in common among these two methods, as well as that of Spouge,
is a free parameter which controls the accuracy of the approximation.
The Lanczos’ method, though, is unique in that the infinite series
term of the formula converges, and the resulting formula defines the
gamma function in right-half planes Re(z) ≥ −r where r ≥ 0 is the
aforementioned free parameter. This is in contrast to the divergent
nature of the series in Stirling’s method.
The main results of this work are

(i) Improved versions of Lanczos’ formulas for computing Γ(z + 1) on

Re(z) ≥ 0. In the specific case of a uniformly bounded relative
error of 10−32 in the right-half plane, a formula which is both
simpler and more efficient than Stirling’s series is found. A one
term approximation similar to Stirling’s formula is also given,
but with a uniform error bound of 0.006 in the right-half plane.
These results stem from the examination of the relative error as
an analytic function of z, and in particular, the behaviour of this
error as |z| → ∞.

(ii) A careful examination of Lanczos’ paper [14], including a proof

of the limit formula stated without proof at the end of the paper.

1.1 Lanczos and His Formula

The object of interest in this work is the formula
√
Γ(z + 1) = 2π (z + r + 1/2)z+1/2 e−(z+r+1/2) Sr (z) (1.1)

where

1 z z(z − 1)
Sr (z) = a0 (r) + a1 (r) + a2 (r) +··· . (1.2)
2 z+1 (z + 1)(z + 2)

This unusual formula is due to Cornelius Lanczos, who published it in

his 11 page 1964 paper A Precision Approximation of the Gamma Func-
tion [14]. The method has been popularized somewhat by its mention
in Numerical Recipes in C [21], though this reference gives no indication

2
Chapter 1. Introduction

as to why it should be true, nor how one should go about selecting trun-
cation orders of the series Sr (z) or values of the parameter r 1 . Despite
this mention, few of the interesting properties of the method have been
explored. For example, unlike divergent asymptotic formulas such as
Stirling’s series, the series Sr (z) converges. Yet, in a manner similar to
a divergent series, the coefficients ak (r) initially decrease quite rapidly,
followed by slower decay as k increases. Furthermore, with increasing
r, the region of convergence of the series extends further and further to
the left of the complex plane, including z on the negative real axis and
not an integer.

Figure 1.1: Cornelius Lanczos

The coefficients ak (r) in (1.2) are expressible in closed form as

functions of the free parameter r ≥ 0, and equation (1.1) is valid for
Re(z + r) > 0, z not a negative integer. Equation (1.1) actually rep-
1 Infact, according to [21], r should be chosen to be an integer, which need
not be the case.

3
Chapter 1. Introduction

resents an infinite family of formulas, one for each value of r, and the
choice of this parameter is quite crucial in determining the accuracy
of the formula when the series Sr (z) is truncated at a finite number
of terms. It is precisely the peculiar r dependency of the coefficients
ak (r), and hence Sr (z) itself, which makes Lanczos’ method interesting.
To get a sense of the r parameter’s influence on Sr (z) at this early
stage, consider the ratios |ak+1 (r)/ak (r)| for sample values r = 1, 4, 7
in Table 1.1. Observe how for r = 1, the relative size of successive

k |ak+1 (1)/ak (1)| |ak+1 (4)/ak (4)| |ak+1 (7)/ak (7)|

0 .31554 1.10590 1.39920
1 .00229 .15353 .33574
2 .32102 .02842 .15091
3 .34725 .00085 .05917
4 .43102 .00074 .01752
5 .50107 .05011 .00258
6 .55757 .36295 .00002
7 .60343 .24357 .00335
8 .64113 .25464 .06382
9 .67258 .27771 .03425
10 .69914 .30151 .58252

Table 1.1: Relative decrease of ak (r), r = 1, 4, 7

coefficients drops very quickly at k = 1, but then flattens out into a

more regular pattern. Compare this with the r = 4 and r = 7 columns
where the drop occurs later, at k = 4 and k = 6 respectively, but the
size of which is increasingly precipitous.
A graphical illustration makes this behaviour more apparent; re-
fer to Figures 1.2, 1.3, and 1.4 for bar graphs of − log |ak+1(r)/ak (r)|,
r = 1, 4, 7, respectively. Observe the large peak to the left of each cor-
responding to the steep drop-off in the coefficients. Considering that
the scale in these graphs is logarithmic, the behaviour is dramatic. No-
tice also that after the initial large decrease, the ratios do not follow a
smooth pattern as in the r = 1 case, but rather, jump around irregu-
larly.
The ak (r) are actually Fourier coefficients of a certain function, and
are bounded asymptotically by a constant times k −2r . However, the ini-

4
Chapter 1. Introduction

12

− log |ak+1 (1)/ak (1)|

−4
0 k 29

Figure 1.2: Relative decrease of ak (1)

12

− log |ak+1 (4)/ak (4)|

−4
0 k 29

Figure 1.3: Relative decrease of ak (4)

tial decrease of the coefficients appears more exponential in nature, and

then switches abruptly to behave according to the asymptotic bound.

5
Chapter 1. Introduction

12

− log |ak+1 (7)/ak (7)|

−4
0 k 29

Figure 1.4: Relative decrease of ak (7)

This peculiar phenomenon is nicely illustrated in Figures 1.5, 1.6 and

1.7 where − log |ak (r)/er | is plotted for k = 0, . . . , 50 using r = 20, 40
and 60, respectively. The abrupt transition in decay (around k ≈ r in
each case) appears unique to Lanczos’ formula, a phenomenon termed
the “Lanczos shelf” in this work. When using the Lanczos formula in
practice, the shelf behaviour suggests that the series (1.2) should be
terminated at about the k ≈ r term, or equivalently, a k-term approx-
imation should use r ≈ k, since the decay rate of coefficients after the
cutoff point slows considerably.
Apart from the observations noted so far, the methods used both
in the derivation of the main formula and for the calculation of the
coefficients makes Lanczos’ work worthy of further consideration. The
techniques used extend to other functions defined in terms of a partic-
ular integral transform.
Lanczos was primarily a physicist, although his contributions to
mathematics and in particular numerical analysis were vast [15]. There
is no question that he had an appreciation for rigour and pure anal-
ysis, and was no doubt adept and skilled technically. At the same
time, however, he seemed to have a certain philosophical view that the
worth of a theoretical result should be measured in part by its practi-

6
Chapter 1. Introduction

300

− log |ak (20)/e20 |

0
0 k 100

Figure 1.5: Lanczos shelf, r = 20

300

− log |ak (40)/e40 |

0
0 k 100

Figure 1.6: Lanczos shelf, r = 40

cality. It is this practicality which is in evidence more so than rigour

in his 1938 paper [11] and his book [12], in which much emphasis is

7
Chapter 1. Introduction

300

− log |ak (60)/e60 |

0
0 k 100

Figure 1.7: Lanczos shelf, r = 60

placed on the benefits of interpolation over extrapolation of analytical

functions (that is, Fourier series versus Taylor series). In these works
one does not find the standard format of theorems followed by proofs
common today. Rather, Lanczos explains his methods in very readable
(and mostly convincing) prose followed by examples, and leaves many
technical details to the reader.
This same style pervades [14], the focus of this study, and is one
of the motivations for the more thorough analysis undertaken here.
A number of statements are given without proof, or even hint of why
they should be in fact be true, and a number of authors have noted this
scarcity of detail. Aspects of his work are variously described as “none
is quite as neat . . . seemingly plucked from thin air” [21], “exceedingly
curious” [29], and “complicated” [27], so a more detailed examination
seems warranted.

1.2 Thesis Overview

The goals of this work are threefold:
(i) To explain in detail Lanczos’ paper, including proofs of all state-

8
Chapter 1. Introduction

ments.

(ii) From a practical point of view, to determine how (1.1) is best

used to compute Γ(z + 1), which reduces to an examination of the
error resulting from truncation of the series (1.2), and to compare
Lanczos’ method with other known methods.

(iii) To note consequences and extensions of Lanczos’ method.

The thesis is organized into eleven chapters, which fall into six more
or less distinct categories. These are

I. History and Background of Γ(z + 1): Chapter 2

II. Lanczos’ Paper: Details and Proofs; Lanczos Limit Formula:

Chapters 3-7

III. Error Discussion: Chapter 8

IV. Comparison of Calculation Methods: Chapter 9

V. Consequences and Extensions of the Theory: Chapter 10

VI. Conclusion and Future Work: Chapter 11

I have attempted to summarize the content and state main results at

the beginning of each chapter, followed by details and proofs in the
body.
The remainder of this introductory chapter is dedicated to an overview
of the six categories noted, with the aim of providing the reader with
a survey of the main ideas, thus serving as a guide for navigating the
rest of the thesis.

History and Background of Γ(z + 1)

The gamma function is a familiar mathematical object for many, being

the standard generalization of the factorial n! for non-negative integer
n. The history of the function is not so well known, however, and so a
brief one is given here, culminating with Euler’s representation
Z ∞
Γ(x) = tx−1 e−t dt
0

9
Chapter 1. Introduction

which is the form with which Lanczos begins his study.

Following the formal definition of Γ(z + 1), a number of standard
identities are given with some explanation but without formal proof.
The chapter concludes with a discussion of computational methods,
beginning with Stirling’s series, which is a generalization of the familiar
Stirling’s formula
√
n! ∼ 2πn nn e−n as n → ∞ .

Stirling’s asymptotic series forms the basis for most computational al-
gorithms for the gamma function and much has been written on its
implementation. An example on the use of Stirling’s series is given
with a short discussion of error bounds.
The second computational method noted is the relatively recent
formula of Spouge [27]. This method is similar in form to Lanczos’
formula but differs greatly in its origin. The formula takes the form
N
" #
X c k
Γ(z + 1) = (z + a)z+1/2 e−(z+a) (2π)1/2 c0 + + ǫ(z)
k=1
z + k

where N = ⌈a⌉ − 1, c0 = 1 and ck is the residue of Γ(z + 1)(z +

a)−(z+1/2) ez+a (2π)−1/2 at z = −k. Although not quite as accurate as
Lanczos’ method using the same number of terms of the series, Spouge’s
method has simpler error estimates, and the coefficients ck are much
easier to compute. This method also extends to the digamma function
Ψ(z + 1) = d/dz [log Γ(z + 1)] and trigamma function Ψ′ (z).

Lanczos’ Paper: Details and Proofs

Chapter 3 examines in detail Lanczos’ paper, beginning with the deriva-

tion of the main formula (1.1). The derivation makes use of Fourier se-
ries in a novel way by first defining an implicit function whose smooth-
ness can be controlled with a free parameter r ≥ 0, and then replacing
that function with its Fourier series. Specifically, beginning with
Z ∞
Γ(z + 1) = tz e−t dt ,
0

10
Chapter 1. Introduction

these steps yield

√
Γ(z + 1/2) = (z + r + 1/2)z+1/2 e−(z+r+1/2) 2

π/2
"√ # (1.3)
2 v(θ)r sin θ
Z
× cos2z θ dθ .
−π/2 log v(θ)

The function v(θ) in the integrand is defined implicitly by

v(1 − log v) = cos2 θ
where v is increasing, and v = 0, 1, e corresponds to θ = −π/2, 0, π/2,
respectively.
Denoting by fr (θ) the term in square brackets in the integrand
of (1.3), and fE,r (θ) its even part, the properties of this last function
are such that it can be replaced by a uniformly convergent Fourier series
∞
a0 (r) X
fE,r (θ) = + ak (r) cos (2kθ) .
2 k=1

Integrating this series against cos2z θ term by term then gives rise to
the series Sr (z) of (1.2) with the help of a handy trigonometric integral.
The derivation of (1.1) can also be carried out using Chebyshev se-
ries instead of Fourier series (really one and the same), which seems
closer in spirit to methods seen later for computing the coefficients
ak (r). In fact, both the Fourier and Chebyshev methods can be gen-
eralized to a simple inner product using Hilbert space techniques, and
this is noted as well.
From the theory of Fourier series, the smoothness of fE,r (θ) governs
the rate of asymptotic decrease of the coefficients ak (r), which in turn
determines the region of convergence of the expansion. This funda-
mental relationship motivates a detailed examination of the properties
of v(θ), fr (θ) and fE,r (θ) in Chapter 4. Practical formulas for these
functions are found in terms of Lambert W functions, and a few rep-
resentative graphs are plotted for several values of r.
Chapter 5 addresses the convergence of the series Sr (z). Once the
growth order of the ak (r) is established, a bound on the functions
(
1 if k = 0,
Hk (z) = z···(z−k+1) (1.4)
(z+1)···(z+k)
if k ≥ 1.

11
Chapter 1. Introduction

appearing in Sr (z) is found which leads to the conclusion that Sr (z)

converges absolutely and uniformly on compact subsets of Re(z) > −r
away from the negative integers, and thus equation (1.1) defines the
gamma function in this region.
Finally, Chapter 6 addresses the practical problem of computing the
ak (r). The standard method provided by Fourier theory, namely the
direct integration
2 π/2
Z
ak (r) = fE,r (θ) cos (2kθ) dθ (1.5)
π −π/2
proves impractical due to the complicated nature of fE,r (θ). Lanczos
overcomes this difficulty with the clever use of Chebyshev polynomials
as follows: define
Fr (z) = 2−1/2 Γ(z + 1/2)(z + r + 1/2)−z−1/2 exp (z + r + 1/2) ,
which is just the integral in (1.3), and recall the defining property of
the 2k th Chebyshev polynomial:
k
X
T2k (cos θ) = C2j,2k cos2j θ
j=0

= cos (2kθ) .
Then (1.5) becomes
π/2 k
2
Z X
ak (r) = fE,r (θ) C2j,2k cos2j θ dθ
π −π/2 j=0

k
2X
= C2j,2k Fr (j) ,
π j=0

a simple linear combination of Fr evaluated at integer arguments.

Aside from Lanczos’ method for computing the coefficients, several
others are noted. In particular, a recursive method similar to Horner’s
rule is described, as well as a matrix method which reduces floating
point operations thus avoiding some round off error. In addition, these
matrix methods provide formulas for the partial fraction decomposition
of the series Sr (z) once terminated at a finite number of terms.

12
Chapter 1. Introduction

The Lanczos Limit Formula

Lanczos concludes his paper with the following formula which he claims
is the limiting behaviour of (1.1), and which defines the gamma function
in the entire complex plane away from the negative integers:
" ∞
#
1 X 2
Γ(z + 1) = 2 lim r z + (−1)k e−k /r Hk (z) . (1.6)
r→∞ 2 k=1

The Hk (z) functions here are as defined by (1.4). He gives no proof, nor
any explanation of why this should be true. In Chapter 7 the limit (1.6)
is proved in detail, a result termed the “Lanczos Limit Formula” in this
work.
The proof is motivated by the plots in Figure 1.8 of the function
fr (θ) from equation (1.3). Notice that as r increases, the area un-
√ 2
2e r=2

fr (θ)

√
2e r=1

√
2 r=0

− π2 0
π
2

Figure 1.8: fr (θ), −π/2 ≤ θ ≤ π/2

13
Chapter 1. Introduction

der fr (θ) appears

√ to be increasingly concentrated near θ = π/2 where
fr (π/2) = 2er .
p
The limit formula (1.6) suggests a rescaling of fr (θ) to r/(2π)fer (θ)e−er ,
which decays rapidly
p to zero on [−π/2, π/2] except at θ = π/2 where
it has the value r/π. p This behaviour is reminiscent of that of the
2
Gaussian distribution r/π exp [−r(θ − π/2) ] on (−∞, π/2], and in-
deed, plots of these two functions in the vicinity of θ = π/2 show a
remarkable fit, even for small values of r. With some work, the limiting
value of (1.3) is shown to be expressible in the form
Z ∞
z
√ 2
Γ(z + 1/2) = lim r cos2z θ re−r(θ−π/2) dθ
r→∞ −∞

which, upon replacement of cos2z θ by its Fourier series, gives rise to

the coefficients via
Z ∞
1 2 2
√ cos (2kθ)r 1/2 e−r(θ−π/2) dt = (−1)k e−k /r .
π −∞

A comparison of (1.6) with the main formula (1.1) suggests that for
large r, r
πr ak (er) 2

er
∼ (−1)k e−k /r , (1.7)
2 e
and plots of these coefficients for various r values supports this rela-
tionship. From the smoothness analysis of fr (θ) it follows that ak (r) =
O k −⌈2r⌉ . Yet, (1.7) suggests a much more rapid decrease of the coef-
ficients for large fixed r and increasing k. Indeed, in his paper, Lanczos
states

Generally, the higher r becomes, the smaller will be the

value of the coefficients at which the convergence begins to
slow down. At the same time, however, we have to wait
longer, before the asymptotic stage is reached.

Although he does not make precise quantitatively how this conclusion is

reached, it seems likely based on the comparison (1.7). Note that this
description is precisely the behaviour of the coefficients observed at
the beginning of this chapter. Unfortunately, the proof of the Lanczos
Limit Formula does not clarify quantitatively Lanczos’ observation, but
it does shed some light on what perhaps motivated him to make it in
the first place.

14
Chapter 1. Introduction

Error Discussion

To use Lanczos’ formula in practice, the series Sr (z) is terminated at

a finite number of terms, and an estimate of the resulting (relative)
error is therefore required. Chapter 8 begins with a survey of existing
error considerations found in the literature, followed by the two prin-
cipal results of the chapter. The first is a simple empirical method for
bounding the relative error. The second is an observation about the
behaviour of the relative error as |z| → ∞. This observation leads to
improved error bounds for certain choices of r as a function of the series
truncation order n.
For the purposes of error estimation, one need only be concerned
with Re(z) ≥ 0 since the reflection formula
πz
Γ(1 + z)Γ(1 − z) =
sin πz
provides estimates for gamma in the left hand plane if the function in
right half plane is known.
For notation, write
Sr (z) = Sr,n (z) + ǫr,n (z)
where Sr,n (z) is the sum of the first n terms and
∞
X
ǫr,n (z) = ak (r)Hk (z) .
k=n+1

Following Lanczos, the function ǫr,n (z) is considered the relative error,
which differs slightly from the standard definition, but which matters
little when it comes to bounding the error. Lanczos gives, for Re(z) ≥ 0,
uniform bounds on |ǫr,n (z)| for seven combinations of n and r values
but does not provide precise details as to how he obtains these figures.
Attempts to reproduce his estimates based on his description were not
successful, and no examination of his estimates appears elsewhere in
the literature. Lanczos’ figures do appear correct, however, based on
the new estimation method which follows.
In order to bound |ǫr,n (z)|, observe that ǫr,n (z) is an analytic func-
tion on Re(z) ≥ 0. It turns out that
n
a0 (r) X
lim ǫr,n (z) = 1 − − ak (r) ,
|z|→∞ 2
k=1

15
Chapter 1. Introduction

a constant denoted ǫ∞ r,n , the error at infinity. Under these conditions,

as a consequence of the maximum modulus principle, the maximum of
|ǫr,n (z)| on Re(z) ≥ 0 must occur on the line Re(z) = 0, possibly at
infinity. Furthermore, thanks to the Schwartz reflection principle, this
maximum will occur on the positive imaginary axis. Thus to bound
|ǫr,n (z)| on Re(z) ≥ 0, it is sufficient to bound |ǫr,n (it)|, 0 ≤ t < ∞:
still a daunting task. The problem is overcome by first transforming the
domain to a finite interval via the mapping z(t) = it/(1 − t), 0 ≤ t < 1,
and then observing that under this mapping the functions Hk (z(t))
have the simple well behaved form
k−1
Y t(i + j) − j
Hk (z(t)) = .
j=0
t(i − j − 1) + j + 1

Thus the maximum of |ǫr,n (z)| on Re(z) ≥ 0 can be easily estimated

empirically by examining the first few terms of |ǫr,n (z(t))| on 0 ≤ t < 1.
This estimation method is used to examine the dependence of ǫr,n (z)
on both the truncation order n and the parameter r. The prevailing
thought in the literature is that n should be chosen as a function of
r; the approach to the problem here is just the opposite: r is selected
as a function of n, with surprising results. Using the relative error in
Stirling’s formula as a guide, the question is: for a given value of n,
can r be chosen so that the relative error at infinity is zero? That is,
does ǫ∞r,n = 0 have solutions? This question motivated an extensive
numerical examination of the relative error functions ǫr,n (z) and ǫ∞ r,n
for n = 0, . . . , 60. The results of this investigation are tabulated in
Appendix C.
Experimentally, ǫ∞r,n = 0 was found to have a finite number of real
solutions. For the values 0 ≤ n ≤ 60 examined, ǫ∞ r,n had at most 2n + 2
real zeros located between r = −1/2 and r = n + 4. Furthermore,
for each n, the largest zero of ǫ∞
r,n was found to give the least uniform
bound on |ǫr,n (z)|. For example, Lanczos gives the following uniform
error bounds for various values of n and r (for Re(z) ≥ 0):

16
Chapter 1. Introduction

n r |ǫr,n (z)| <

1 1 0.001
1 1.5 0.00024
2 2 5.1 × 10−5
3 2 1.5 × 10−6
3 3 1.4 × 10−6
4 4 5 × 10−8
6 5 2 × 10−10

By contrast, selecting r to be the largest zero of ǫ∞

r,n yields the following
dramatically improved uniform error bounds (again, for Re(z) ≥ 0):

n r |ǫr,n (z)| <

1 1.489194 1.0 × 10−4
2 2.603209 6.3 × 10−7
3 3.655180 8.5 × 10−8
4 4.340882 4.3 × 10−9
6 6.779506 2.7 × 10−12

Comparison of Calculation Methods

The purpose of this chapter is to compare the methods of Lanczos,

Spouge and Stirling in the context of an extended precision compu-
tation. For each method, a detailed calculation of Γ(20 + 17i) with
relative error |ǫρ | < 10−32 is carried out. Additionally, formulas for
Γ(1 + z) with uniform error bound of 10−32 are given for each.
The conclusion is that each method has its merits and shortcomings,
and the question of which is best has no clear answer. For a uniformly
bounded relative error, Lanczos’ method seems most efficient, while
Stirling’s series yields very accurate results for z of large modulus due
to its error term which decreases rapidly with increasing |z|. Spouge’s
method, on the other hand, is the easiest to implement thanks to its
simple formulas for both the series coefficients and the error bound.

Consequences and Extensions of the Theory

Chapter 10 discusses a variety of results which follow from various as-

pects of Lanczos’ paper and the theory of earlier chapters.

17
Chapter 1. Introduction

The main result is a generalization of the integral transform in (1.3)

which is termed a “Lanczos transform” in this work. Briefly, if F (z) is
defined for Re(z) ≥ 0 as
Z π/2
F (z) = cos2z θg(θ) dθ
−π/2

where g(θ) ∈ L2 [−π/2, π/2] is even, then

√ Γ(z + 1/2) 1

z z(z − 1)
F (z) = π a0 + a1 + a2 +··· .
Γ(z + 1) 2 z+1 (z + 1)(z + 2)

The ak in the series are the Fourier coefficients of g(θ), and these are
given by
k
2X
ak = C2j,2k F (j) ,
π j=0

exactly as in the gamma function case. Again, just as in the gamma

function case, the smoothness of g has a direct influence on the domain
and speed of convergence of the series. The relative error resulting from
truncating the series at a finite number of terms is constant at infin-
ity, and so the empirical error bounding methods used in the gamma
function case also apply to the more general Lanczos transform.
Aside from the Lanczos transform, two non-trivial combinatorial
identities which follow directly from the Lanczos Limit Formula are
noted. The third result is a variation on Stirling’s formula. A result of
the work in previous chapters was the one term approximation
√
Γ(z + 1) = (z + r + 1/2)(z+1/2) e−(z+r+1/2) 2π (1 + ǫr,0 (z)) ,
.
where r = 0.319264, the largest zero of ǫ∞
r,0 , and |ǫr,0 (z)| < 0.006 every-
where in the right-half plane Re(z) ≥ 0. The natural question is: can
r be chosen as a function of z so that ǫr(z),0 (z) = 0, thus yielding
√
Γ(z + 1) = (z + r(z) + 1/2)(z+1/2) e−(z+r(z)+1/2) 2π ?

The answer is shown to be yes, with the help of Lambert W functions,

and numerical checks indicate that the functions r(z) vary little over
the entire real line.

18
Chapter 1. Introduction

Future Work and Outstanding Questions

The concluding chapter notes a number of outstanding questions and

areas requiring further investigation. This topic is broken down into
three categories: unresolved problems from Lanczos’ paper itself, ques-
tions arising out of the theory developed in this study, particularly
with respect to error bounds, and finally, a conjectured deterministic
algorithm for calculating the gamma function based on the numerical
evidence of Appendix C.
The last of these is of greatest practical interest. Letting r(n) denote
the largest zero of ǫ∞ r,n and Mr(n),n the maximum of |ǫr,n (it)| for 0 ≤
t < ∞, a plot of (n, − log Mr(n),n ) reveals a near perfect linear trend
n = −a log Mr(n),n + b. The algorithm is then clear: given ǫ > 0, select
n = ⌈−a log ǫ + b⌉, and set r = r(n). Then (1.1) truncated after the
nth term computes Γ(1 + z) with a uniformly bounded relative error
of at most ǫ. Although a proof is not available at present, the data is
compelling and further investigation is warranted.

19
Chapter 2
A Primer on the Gamma
Function

This chapter gives a brief overview of the history of the gamma func-
tion, followed by the definition which serves as the starting point of
Lanczos’ work. A summary of well known results and identities involv-
ing the gamma function is then given, concluding with a survey of two
computational methods with examples: Stirling’s series and Spouge’s
method.
The introductory material dealing with the definition and resulting
identities is standard and may be safely skipped by readers already
familiar with the gamma function. The computational aspects of the
function, on the other hand, may be less familiar and readers may
therefore wish to include the material starting with Section 2.5.

2.1 First Encounters with Gamma

Most any student who has taken a mathematics course at the senior
secondary level has encountered the gamma function in one form or
another. For some, the initial exposure is accidental: n! is that button
on their scientific calculator which causes overflow errors for all but
the first few dozen integers. The first formal treatment, however, is
typically in the study of permutations where for integer n ≥ 0 the
factorial function first makes an appearance in the form
(
1 if n = 0,
n! =
n(n − 1)! else.

20
Chapter 2. A Primer on the Gamma Function

Once some familiarity with this new object is acquired (and students
are convinced that 0! should indeed equal 1), it is generally applied to
counting arguments in combinatorics and probability. This is the end
of the line for some students. Others continue on to study calculus,
where n! is again encountered in the freshman year, typically in the
result
dn n
x = n!
dxn
for n a non-negative integer. Closely related to this result is the ap-
pearance of n! in Taylor’s formula. Some are asked to prove
Z ∞
n! = tn e−t dt (2.1)
0

as homework, while a fortunate few learn how (2.1) “makes sense” for
non-integer n. Some even go so far as to prove
√
n! ∼ 2πn nn e−n as n → ∞ .

At this point the audience begins to fracture: some leave math-

ematics and the gamma function forever, having satisfied their credit
requirement. Others push on into engineering and the physical sciences,
and a few venture into the more abstract territory of pure mathemat-
ics. The last two of these groups have not seen the end of the gamma
function. For the physical scientists, the next encounter is in the study
of differential equations and Laplace transform theory, and for those
who delve into asymptotics, computation of the gamma function is a
classic application of the saddle point method. In statistics the gamma
function forms the basis of a family of density functions which includes
the exponential and chi-square distributions.
For the mathematicians the show is just beginning. In complex
analysis they learn that not only can the gamma function be extended
to non-integer arguments, but even to an analytic function away from
the negative integers. Computationally, they learn about Stirling’s for-
mula and the associated Stirling series. In number theory, the gamma
function is the springboard for developing the analytic continuation
of Riemann’s zeta function ζ(s), which leads into the prime number
theorem and a tantalizing glimpse of the Riemann Hypothesis.
The gamma function has many guises, but exactly how is it defined,
and more importantly for our purposes, how does one compute it? So

21
Chapter 2. A Primer on the Gamma Function

far, a precise definition of the function, and notation for it (save the
n! used) has been avoided. The gamma function is a generalization of
n!, but the question remains: what definition generalizes the factorial
function in the most natural way?

2.2 A Brief History of Gamma

Most works dealing with the gamma function begin with a statement
of the definition, usually in terms of an integral or a limit. These def-
initions generalize the factorial n!, but there are many other functions
which interpolate n! between the non-negative integers; why are the
standard definitions the “right” ones?
The life of gamma is retraced in the very entertaining paper of
Davis [6], according to which the year 1729 saw “the birth” of the
gamma function as Euler studied the pattern 1, 1 · 2, 1 · 2 · 3, . . . .
The problem was simple enough: it was well known that interpolating
formulas of the form
n(n + 1)
1+2+···+n =
2
existed for sums; was there a similar formula f (n) = 1 · 2 · · · n for prod-
ucts? The answer, Euler showed, was no; the products 1 · 2 · · · n would
require their own symbol. The notation n! was eventually adopted
(though not universally) to denote the product, however this notation
did not arrive on the scene until 1808 when C. Kramp used it in a
paper. Today the symbol n! is normally read “n-factorial”, although
some early English texts suggested the somewhat amusing reading of
“n-admiration”. Refer to [4] for a detailed history of the notation.
Euler went on further and found representations of the factorial
which extended its domain to non-integer arguments. One was a prod-
uct
m!(m + 1)x
x! = lim , (2.2)
m→∞ (x + 1) · · · (x + m)

and the second an integral representation:

Z 1
x! = (− log t)x dt . (2.3)
0

22
Chapter 2. A Primer on the Gamma Function

The product (2.2) converges for any x not a negative integer. The
integral (2.3) converges for real x > −1. The substitution u = − log t
puts (2.3) into the more familiar form
Z ∞
x! = ux e−u du .
0

In 1808 Legendre introduced the notation which would give the

gamma function its name1 . He defined
Z ∞
Γ(x) = tx−1 e−t dt , (2.4)
0

so that Γ(n + 1) = n!. It is unclear why Legendre chose to shift the

argument in his notation. Lanczos remarks in the opening paragraph
of his paper [14]:

The normalization of the gamma function Γ(n + 1) in-

stead of Γ(n) is due to Legendre and void of any rationality.
This unfortunate circumstance compels us to utilize the no-
tation z! instead of Γ(z + 1), although this notation is ob-
viously highly unsatisfactory from the operational point of
view.

Lanczos’ was not the only voice of objection. Edwards, in his well-
known treatise on the Riemann zeta function [8], reverts to Gauss’
notation Π(s) = s!, stating:

Unfortunately, Legendre subsequently introduced the no-

tation Γ(s) for Π(s − 1). Legendre’s reasons for considering
(n−1)! instead of n! are obscure (perhaps he felt it was more
natural to have the first pole occur at s = 0 rather than at
s = −1) but, whatever the reason, this notation prevailed
in France and, by the end of the nineteenth century, in the
rest of the world as well. Gauss’s original notation appears
to me to be much more natural and Riemann’s use of it
gives me a welcome opportunity to reintroduce it.

It is true, however, that (2.4) represents Γ(x) as the Mellin transform

of exp (−t), so that Legendre’s normalization is the right one in some
1 Whittaker and Watson [30] give 1814 as the date, however.

23
Chapter 2. A Primer on the Gamma Function

sense, though Mellin was not born until some twenty one years after
Legendre’s death. An interesting recent investigation into the history
of the Γ notation appears in the paper by Gronau [10].
The question remains, are Euler’s representations (2.2) and (2.3)
the “natural” generalizations of the factorial? These expressions are
equivalent for x > −1 (see [16]), so consider the integral form (2.3),
or equivalently, (2.4). Some compelling evidence in support of (2.4)
comes in the form of the complete characterization of Γ(x) by the Bohr-
Mollerup theorem (see [23]):
Theorem 2.1. Suppose f : (0, ∞) → (0, ∞) is such that f (1) = 1,
f (x + 1) = xf (x), and f is log-convex. Then f = Γ.

So as a function of a real variable, Γ(x) has nice behaviour. But if

gamma is extended beyond the real line, its properties are R ∞even more
z−1 −t
satisfying: as a function of a complex variable z, Γ(z) = 0 t e dt
defines an analytic function on Re(z) > 0. Along with the fundamental
recursion Γ(z + 1) = zΓ(z), Γ can then be extended to an analytic
function on z ∈ C \ {0, −1, −2, −3, . . .}. Indeed, Euler’s formula (2.2)
converges and defines Γ on this same domain.
The definition created to describe a mathematical object is an arti-
ficial absolute which, it is hoped, precisely captures and describes every
aspect which characterizes the object. However, it is impossible to state
with certainty that a particular definition is in some way the correct
one. In the present case, a definition of the correct generalization of
n! is sought. The Bohr-Mollerup theorem and the analytic properties
of Euler’s Γ strongly suggest that (2.4) is indeed the right definition.
The appearance of Euler’s Γ in so many areas of mathematics, and its
natural relationship to so many other functions further support this
claim. This viewpoint is perhaps best summarized by Temme in [28]:
This does not make clear why Euler’s choice of general-
ization is the best one. But, afterwards, this became abun-
dantly clear. Time and again, the Euler gamma function
shows up in a very natural way in all kinds of problems.
Moreover, the function has a number of interesting proper-
ties.
R∞
So is Γ(z) = 0 tz−1 e−t dt the natural extension of n! ? It seems so,
and we will take it as such and study its properties, and in particular,
how to compute it.

24
Chapter 2. A Primer on the Gamma Function

2.3 Definition of Γ
As noted, there are several equivalent ways to define the gamma func-
tion. Among these, the one which served as the starting point for
Lanczos’ paper [14] will be used:
Definition 2.1. For z ∈ C with Re(z) > −1, the gamma function is
defined by the integral
Z ∞
Γ(z + 1) = tz e−t dt . (2.5)
0

From this definition, we immediately state

Theorem 2.2. For Γ(z + 1) defined by (2.5):

(i) Γ(z + 1) is analytic on Re(z) > −1;

(ii) Γ(z + 1) = zΓ(z) on Re(z) > 0 ;

(iii) Γ(z + 1) extends to an analytic function on C with simple poles

at the negative integers.

Proof of Theorem 2.2:

R∞
(i) Let Ω denote the domain Re(z) > −1. The integrand of 0 tz e−t dt
is continuous as a function of t and z, and for fixed t > 0 is ana-
lytic on Ω. Further the integral converges uniformly on compact
subsets of Ω. Hence Γ(z + 1) is analytic on Ω.

(ii) For Re(z) > 0, integrating by parts gives

Z ∞
Γ(z + 1) = tz d(−e−t )
0
Z ∞
z −t ∞
= −t e 0 + z tz−1 e−t dt
0
= zΓ(z)

(iii) For Re(z) > 0 and any integer n ≥ 0, the fundamental recursion
Γ(z + 1) = zΓ(z) implies
Γ(z + n + 1)
Γ(z) = Qn . (2.6)
k=0 (z + k)

25
Chapter 2. A Primer on the Gamma Function

The right hand side of (2.6) is analytic on Re(z) > −n − 1, z not

a negative integer, and is equal to Γ(z) on Re(z) > 0. Hence (2.6)
is the unique analytic continuation of Γ(z) to the right half plane
Re(z) > −n−1, z not a negative integer. Since n was an arbitrary
non-negative integer, Γ(z) can be extended to an analytic function
on C \ {0, −1, −2, −3, . . .}.
For z = −n a negative integer, (2.6) shows that Γ has a simple
pole with residue

Γ(z + n + 1)
lim (z + n)Γ(z) = lim (z + n) Qn
k=0 (z + k)
z→−n z→−n

Γ(1)
= Qn−1
k=0 (−n + k)
1
= Qn−1
k=0 (−n + k)
(−1)n
= .
n!

An interesting generalization of (2.5) due to Cauchy leads to another

analytic continuation of the gamma function. Let σ = Re(−z − 1).
Then Z ∞ !
X (−t)k
Γ(z + 1) = tz e−t − dt . (2.7)
0 0≤k<σ
k!

An empty sum in the integrand of (2.7) is to be considered zero.

This integral converges and is analytic for z ∈ C such that Re(z) 6=
−1, −2, −3, . . . and reduces to (2.5) if Re(z) > −1. See [30, pp.243–
244]
Figure 2.1 shows a graph of Γ(z + 1) for z real. Note the the poles
at z = −1, −2, . . . and the contrast between the distinctly different
behaviours of the function on either side of the origin.
Figure 2.2 shows a surface plot of |Γ(x + iy)| which illustrates the
behaviour of the function away from the real axis.

26
Chapter 2. A Primer on the Gamma Function

Γ(z + 1)
3

−2 2 z

−3

Figure 2.1: Γ(z + 1), −4 < z < 4

2.4 Standard Results and Identities

From the proof of Theorem 2.2 it is clear that the fundamental recursion
Γ(z + 1) = zΓ(z) holds for all z ∈ C \ {0, −1, −2, −3, . . .}. There are
many more properties of the gamma function which follow from the
definition and Theorem 2.2. The main ones are stated here without
proof; the interested reader is referred to [16, pp. 391–412] for details.
The first result is Euler’s formula
n!nz
Γ(z) = lim (2.8)
n→∞ z(z + 1) · · · (z + n)

which converges for all z ∈ C\{0,

R n −1, −2,
n −3, . . .}. This identity follows
from integrating by parts in 0 1 − nt tz−1 dt to get
Z n n
t n!nz
1− tz−1 dt =
0 n z(z + 1) · · · (z + n)

27
Chapter 2. A Primer on the Gamma Function

0
2
−2
0
0
2 y
x −2

Figure 2.2: |Γ(x + iy)|, −3 < x < 4, −3 < y < 3

and then letting n → ∞ .

Rewriting (2.8) in reciprocal form as
n
z(z + 1) · · · (z + n) z(− log n+ n
P
k=1 1/k)
Y z −z/k
lim = lim ze 1 + e
n→∞ n!nz n→∞
k=1
k

yields the Weierstrass product

∞
1 γz
Y z −z/k
= ze 1+ e . (2.9)
Γ(z) k=1
k

The γ appearing in (2.9)

Pnis the Euler (or sometimes Euler- Mascheroni)
constant γ = limn→∞ k=1 (1/k) − log n ≈ 0.5772 · · · .
Writing both Γ(z) and Γ(−z) using (2.9) and recalling the infinite

28
Chapter 2. A Primer on the Gamma Function

product representation of the sine function

∞
z2
Y
sin z = z 1− 2 2 ,
k=1
k π

we arrive at the very useful reflection formula

πz
Γ(1 + z)Γ(1 − z) = . (2.10)
sin πz
Equation (2.10) is of great practical importance since it reduces the
problem of computing the gamma function at arbitrary z to that of z
with positive real part.

2.5 Stirling’s Series and Formula

Among the well known preliminary results on the gamma function, one
is of particular importance computationally and requires special atten-
tion: Stirling’s series. This is an asymptotic expansion for Γ(z + 1)
which permits evaluation of the function to any prescribed accuracy,
and variations of this method form the basis of many calculation rou-
tines.

2.5.1 The Formula

Stirling’s series is named after James Stirling (1692–1770) who in 1730
published the simpler version for integer arguments commonly known
as Stirling’s Formula. The more general Stirling’s series can be stated
thus:
" N
#
Y
log Γ(z + 1) (z + k) = (z + N + 1/2) log (z + N)
k=1
1
− (z + N) + log 2π
2
n (2.11)
X B2j
+
j=1
2j(2j − 1)(z + N)2j−1
Z ∞
B2n (x)
− dx .
0 2n(z + N + x)2n

29
Chapter 2. A Primer on the Gamma Function

The product in the left hand side of (2.11) is understood to be one

if N = 0. The B2j are the Bernoulli numbers, which are defined by
writing t/(et − 1) as a Maclaurin series:
∞
t X Bj
t
= tj .
e −1 j=0
j!

The function B2n (x) is the 2nth Bernoulli polynomial defined over R by
periodic extension of its values on [0, 1].
By judiciously selecting n and N the absolute error in log Γ(z + 1),
that is, the absolute value of the integral term in (2.11), can be made as
small as desired, and hence the same can be achieved with the relative
error of Γ(z + 1) itself. Equation (2.11) is valid in the so-called slit
plane C \ {t ∈ R | t ≤ −N}.
Equation (2.11) with N = 0 is derived using Euler-Maclaurin sum-
mation. Alternatively, again with N = 0, begin with the formula of
Binet obtained from an application of Plana’s theorem,
Z ∞
1 arctan (t/z)
log Γ(z + 1) = (z + 1/2) log z − z + log 2π + 2 dt ,
2 0 e2πt − 1
and expand arctan (t/z) in the last term into a Taylor polynomial
with remainder. The resulting series of integrals produce the terms
of (2.11) containing the Bernoulli numbers [30, pp.251–253]. The gen-
eral form (2.11) with N ≥ 1 is obtained by replacing z with z + N and
applying the fundamental recursion.
Exponentiation of Stirling’s series (with N = 0) yields Stirling’s
asymptotic formula for Γ(z + 1) itself:

−z z+1/2 1/2 1 1
Γ(z + 1) = e z (2π) 1+ + +··· . (2.12)
12z 288z 2
Equation (2.12) also follows from applying Laplace’s method (or equiva-
lently, the saddle point method) to (2.5). Truncating the series in (2.12)
after the constant term yields the familiar Stirling’s formula noted in
the introduction:
√
Γ(z + 1) ∼ 2πz z z e−z as |z| → ∞ . (2.13)

It may be tempting to let n → ∞ in (2.11) with the hope that the

integral error term will go to zero, but unfortunately the series contain-
ing the B2j becomes divergent. The terms of this sum initially decrease

30
Chapter 2. A Primer on the Gamma Function

in modulus, but then begin to grow without bound as a result of the

rapid growth of the Bernoulli numbers. The odd Bernoulli numbers are
all zero except for B1 = −1/2, while for the even ones, B0 = 1, and for
j ≥ 1,
(−1)j+1 2(2j)!ζ(2j)
B2j = . (2.14)
(2π)2j
Since ζ(2j) = ∞ −2j
→ 1 rapidly with increasing j, B2j = O ((2j)!(2π)−2j )
P
k=1 k
(see [8, p.105]).
What is true, however, is that the integral error term in (2.11) tends
to zero uniformly as |z| → ∞ in any sector | arg (z + N)| ≤ δ < π. In
terms of Γ(z + 1) itself, this means the relative error tends to zero
uniformly as |z| → ∞ in any sector | arg (z + N)| ≤ δ < π.
In order to use (2.11) to evaluate the gamma function, an estimate
of the error (as a function of n, N, and z) is required when the integral
term of (2.11) is omitted. The following result is due to Stieltjes [8,
p.112]:

Theorem 2.3. Let θ = arg (z + N), where −π < θ < π, and denote
the absolute error by
Z ∞
B2n (x)
EN,n (z) = dx .
0 2n(z + N + x)2n

Then
2n+2
1 B2n+2
|EN,n (z)| ≤ (2n + 2)(2n + 1)(z + N)2n+1 . (2.15)

cos (θ/2)

In other words, if the series in (2.11) is terminated after the B2n

term and the integral term is omitted, the resulting error is at most
(cos (θ/2))−2n−2 times the B2n+2 term.
To get an idea of how to use Stirling’s series and the error esti-
mate (2.15), consider the following

Example 2.1. Compute Γ(7+13i) accurate to within an absolute error

of ǫ = 10−12 .

Solution: For convenience, let EN,n denote the integral error term
EN,n (6 + 13i) in this example. The first task is to translate the absolute

31
Chapter 2. A Primer on the Gamma Function

bound ǫ = 10−12 on the error in computing Γ(7 + 13i) into a bound

on EN,n . Let log (GPN ) denote the left hand side of (2.11), where G
represents the true value of Γ(7 + 13i) and PN the product term. Let
log GN,n denote the right hand side of (2.11) without the integral error
term. Then (2.11) becomes

log (GPN ) = log GN,n + EN,n ,

so that
GN,n EN,n
G= e ,
PN
and the requirement is

GN,n E G N,n
PN e −
N,n <ǫ.
PN

That is, it is desired that

E
e N,n − 1 < ǫ
|GN,n /PN |

from which
2
EN,n

EN,n ǫ
|EN,n | 1 + + + · · · < .
2! 3! |GN,n /PN |

The infinite series is at most e in modulus provided EN,n ≤ 1, so that if

|EN,n | < ǫ/(e|GN,n /PN |), the prescribed accuracy is guaranteed. This
last estimate is self-referencing, but GN,n /PN can be estimated using
Stirling’s formula (2.13). That is, we require
ǫez−1

|EN,n | < √

2πz z+1/2

≈ 5 × 10−12

upon setting z = 6 + 13i.

Now that a target bound for |EN,n | has been determined, n and N
must now be selected to meet the target. In Table 2.1 are listed bounds
on |EN,n | as given by (2.15) for z = 6 + 13i and various combinations
of (N, n).

32
Chapter 2. A Primer on the Gamma Function

n N =0 N =1 N =2 N =3 N =4
1 1.9 × 10−6 1.6 × 10−6 1.3 × 10−6 1.1 × 10−6 9.7 × 10−7
2 3.7 × 10−9 2.8 × 10−9 2.2 × 10−9 1.7 × 10−9 1.3 × 10−9
3 1.9 × 10−11 1.3 × 10−11 9.1 × 10−12 6.4 × 10−12 4.4 × 10−12
4 1.9 × 10−13 1.2 × 10−13 7.3 × 10−14 4.6 × 10−14 2.9 × 10−14
5 2.9 × 10−15 1.6 × 10−15 9.3 × 10−16 5.3 × 10−16 3.1 × 10−16

Table 2.1: Upper bound on |EN,n (6 + 13i)| in Stirling Series

We can see that N = 0 with n = 4 is sufficient to produce the

desired accuracy. Using these values with z = 6 + 13i, and leaving off
the error term in (2.11) then gives

log Γ(7 + 13i) ≈ log G0,4

1
= (6 + 13i + 1/2) log (6 + 13i) − (6 + 13i) + log 2π
2
4
X B2j
+
j=1
2j(2j − 1)(6 + 13i)2j−1
.
= −2.5778902638380984 + i 28.9938056395651838 .

Taking exponentials yields the final approximation

.
Γ(7 + 13i) ≈ G0,4 = −0.0571140842611710 − i 0.0500395762571980 .

Using the symbolic computation package Maple with forty digit preci-
.
sion, the absolute error in the estimate is found to be |Γ(7 + 13i) − G0,4 | =
2.5 × 10−15 , which is well within the prescribed bound.
Observe from Table 2.1 that N = 4 with n = 3 would also have
given the desired accuracy, in which case only n = 3 terms of the series
would be needed.

The reader has no doubt noticed the ad-hoc manner in which n and
N were selected in this example. The question of how best to choose n
and N given z and absolute error ǫ is not an easy one, and will not be
addressed here.

33
Chapter 2. A Primer on the Gamma Function

It is a simpler matter, however, to determine the maximum accuracy

that can be expected for a given choice of N, both for a fixed z and
uniformly, and this is done in the following sections.

2.5.2 Some Remarks about Error Bounds

From Table 2.1 it may seem difficult to believe that as n continues
to increase, the error will reach a minimum and then eventually grow
without bound. However, for a fixed z and N it is possible to determine
the value of n which gives the least possible error bound that can be
achieved with (2.15). More generally, for each N ≥ 1, it is possible to
establish the best possible uniform error estimate reported by (2.15)
for Re(z) ≥ 0.
In the following discussion, let
2n+2
1 B2n+2
UN,n (z) = ,
cos (θ/2) (2n + 2)(2n + 1)(z + N)2n+1
which is the right hand side of (2.15). Recall that θ = arg (z + N).

Minimum of |UN,n(6 + 13i)|

For each N, what is the best accuracy that can be achieved in Exam-
ple 2.1 using the error estimate in Stirling’s Series (2.15)? Throughout
this section let |UN,n | = |UN,n (6 + 13i)|. To find the minimum value of
|UN,n |, with N fixed, it is sufficient to determine the n value at which
the sequence of |UN,n | changes from decreasing to increasing, which cor-
responds to the first n such that |UN,n /UN,n+1 | ≤ 1. Using the explicit
form of the Bernoulli numbers (2.14), the task then is to find the least
n such that
(2π)2 (z + N)2 cos2 (θ/2) ζ(2n + 2)

≤1.
(2n + 2)(2n + 1) ζ(2n + 4)
That is, we want the least n satisfying
s
3 1 ζ(2n + 4)
|π(z + N) cos (θ/2)| ≤ n 1+ + 2 .
2n 2n ζ(2n + 2)
√
The square root term is at most 14π/7 when n = 1 and decreases to
1 rapidly, so that if n is chosen
n ≈ ⌈π |z + N| cos (θ/2)⌉ (2.16)

34
Chapter 2. A Primer on the Gamma Function

a minimal error should result.

q √
In the case of z = 6 + 13i, cos (θ/2) = 1/2 + 3/ 205 and N = 0,
(2.16) predicts a minimum error for n = 38. Using Maple the cor-
responding error is found to be |E0,38 | < 1.3 × 10−34 . This bound is
impressive, but bear in mind that the corresponding Bernoulli number
required to achieve this is |B76 | ≈ 8 × 1050 . This is the best possible
accuracy for z = 6 + 13i with N = 0; greater accuracy is possible as
N is increased from zero. For example, with N = 4, (2.16) predicts a
minimum error for n = 47, at which |E4,47 | < 6.7 × 10−42 .

Uniform Bound of |UN,n(z)|

As demonstrated, for fixed z and N there is a limit to the accuracy

which can be achieved by taking more and more terms of the Stirling
series. It is worthwhile to determine this limiting accuracy in the form
of the best possible uniform error bound as a function of N ≥ 0. By
the reflection formula (2.10), it is enough to consider only the right half
plane Re(z) ≥ 0. Now for N ≥ 0 fixed and z + N = u + iv, |UN,n (z)|
will be worst possible where

| cos (θ/2)|2n+2 |u + iv|2n+1

is a minimum in u ≥ N. If N = 0, letting u + iv → 0 along the positive

real axis shows that |U0,n | grows
p without bound. Assume then √ that
N ≥ 1. Writing cos (θ/2) = (1 + cos θ)/2 and cos θ = u/ u2 + v 2
puts (2.5.2) in the form
1 n/2 √ n+1
| cos (θ/2)|2n+2 |u + iv|2n+1 = u2 + v 2 u2 + v 2 + u
2n+1
which is clearly minimized when u and v are minimized, that is, at
v = 0 and u = N. In terms of z and θ, this says that |UN,n (z)| is worst
possible at z = θ = 0. By (2.16), n should then be chosen n ≈ ⌈πN⌉,

35
Chapter 2. A Primer on the Gamma Function

and we have

B2n+2
|UN,n | ≤

2n+1

(2n + 2)(2n + 1)N

2 (2n + 2)! ζ(2n + 2)

= from (2.14)
(2π)2n+2 (2n + 2) (2n + 1) N 2n+1

2 (2n)!
≈ since ζ(2n + 2) ≈ 1
(2π)2n+2 N 2n+1

√
2 4πn (2n)2n e−2n
≈ from (2.13)
(2π)2n+2 N 2n+1

e−2πN
= √ upon setting n ≈ πN .
π N

So selecting n ≈ ⌈πN⌉ results in a uniform bound which decreases

rapidly with N. The problem remains, however, that the larger N,
and hence n becomes, so too do the Bernoulli numbers required in the
Stirling series. In practice, unless one is dealing with z values near the
origin, one can use much smaller values of N which result in acceptable
error bounds.

2.6 Spouge’s Method

To conclude this survey of standard results, special mention is made of
the 1994 work of Spouge [27]. In that work, the author develops the
formula
" N
#
X c k (a)
Γ(z + 1) = (z + a)z+1/2 e−(z+a) (2π)1/2 c0 + + ǫ(z) (2.17)
k=1
z + k

which is valid for Re(z + a) > 0. The parameter a is real, N = ⌈a⌉ − 1,

c0 (a) = 1 and ck (a) is the residue of Γ(z + 1)(z + a)−(z+1/2) ez+a (2π)−1/2
at z = −k. Explicitly, for 1 ≤ k ≤ N, this is
1 (−1)k−1
ck (a) = √ (−k + a)k−1/2 e−k+a . (2.18)
2π (k − 1)!

36
Chapter 2. A Primer on the Gamma Function

Spouge’s formula has the very simple relative error bound ǫS (a, z)

ǫ(z)
|ǫS (a, z)| =

Γ(z + 1)(z + a)−(z+1/2) ez+a (2π)−1/2

√
a 1
< a+1/2
,
(2π) Re(z + a)

provided a ≥ 3. Thus for z in the right half plane Re(z) ≥ 0, ǫS (a, z)

has the uniform bound
1
|ǫS (a, z)| < √ . (2.19)
a(2π)a+1/2

Though similar in form to Lanczos’ formula (note for example, the

free parameter a), Spouge’s work differs greatly in the derivation, mak-
ing extensive use of complex analysis and residue calculus.
To see how Spouge’s method works in practice, we revisit the com-
putation of Γ(7 + 13i) from Example 2.1:
Example 2.2. Estimate Γ(7+13i) accurate to within an absolute error
of ǫ = 10−12 using Spouge’s method.

Solution: An absolute error bound of ǫ < 10−12 means the relative

error must be bounded by

ǫ
|ǫS (a, z)| <

Γ(7 + 13i)

−12 6+13i

10 e
≈ √ by Stirling’s formula
2π(6 + 13i)6+13i+1/2

.
= 1.3 × 10−11 .

By plotting |ǫS (a, 6 + 13i)|, a = 12.5 and hence N = 12 terms of the

series (2.17) are sufficient to achieve this bound. The calculation of
the coefficients (2.18) for these values yields Table 2.2. Equation (2.17)

37
Chapter 2. A Primer on the Gamma Function

k ck (a)
0 +1.0000000000000 × 100
1 +1.3355050294248 × 105
2 −4.9293093529936 × 105
3 +7.4128747369761 × 105
4 −5.8509737760400 × 105
5 +2.6042527033039 × 105
6 −6.5413353396114 × 104
7 +8.8014596350842 × 103
8 −5.6480502412898 × 102
9 +1.3803798339181 × 101
10 −8.0781761698951 × 10−2
11 +3.4797414457425 × 10−5
12 −5.6892712275042 × 10−12

Table 2.2: Coefficients of Spouge’s Series, a = 12.5

then gives Γ(7 + 13i) ≈ −0.0571140842611682 − i0.0500395762571984,

which differs in absolute value from Γ(7 + 13i) by less than 2.9 × 10−16 .
Compared to Stirling’s series, many more terms of Spouge’s series
are required to achieve the same accuracy, but the individual terms are
much easier to compute, and the selection criterion for a and hence N
is straightforward.

Spouge’s work is important for several reasons. The first is that

the coefficients ck (a) are simpler to compute than those of the Lanczos
series. The second is that Spouge gives a simpler yet more accurate
version of Stirling’s formula. And finally, Spouge’s approximation and
error estimates apply not only to Γ(z + 1), but also to the digamma
function Ψ(z + 1) = d/dz [log Γ(z + 1)] and trigamma function Ψ′ (z).
To see the link between (2.17) and Lanczos’ formula (1.1), write
a = r + 1/2 and resolve the first N terms of the series (1.2) into partial

38
Chapter 2. A Primer on the Gamma Function

fractions,

1 z z · · · (z − N + 1)
a0 (r) + a1 (r) + · · · + aN (r) + ǫ(z)
2 z+1 (z + 1) · · · (z + N)
N
X bk (r)
= b0 (r) + + ǫ(z) .
z+k
k=1

Comparing this with (2.17), the bk (r) obtained from truncating the
Lanczos series are the approximate residues of Γ(z+1)(z+a)−(z+1/2) ez+a (2π)−1/2
at z = −k, and the larger N becomes the better the approximation.

2.7 Additional Remarks

The standard definition and properties of the gamma function can be
found in many works dealing with special functions and analysis. For
this study the especially clear treatments given in [28, pp.41–77] and
[23, pp.192–195] are of note. The references [16, pp. 391–412] and
[30, pp.235–264] are also worthy of mention. A thorough treatment of
Stirling’s series can be found in [8, p.105]. Finally, the work of Artin [2]
is of interest for the treatment of Γ(x) for x real from the point of view
of convexity.
For a survey and comparison of various computational methods for
computing the gamma function, including Stirling’s series, see the paper
by Ng [20]. Among the methods covered there is that of Spira [26] in
which a simpler error bound on Stirling’s series than (2.15) is given.

39
Chapter 3
The Lanczos Formula

We arrive at last at the examination of Lanczos’ paper itself, beginning

in this chapter with the derivation of the main formula (1.1). The
derivation consists of three steps, which in a nutshell can be broken
down as
Z ∞
Γ(z + 1/2) = tz−1/2 e−t dt
0

√
= (z + r + 1/2)z+1/2 e−(z+r+1/2) 2
Z e
× [v(1 − log v)]z−1/2 v r dv (step I)
0

√
= (z + r + 1/2)z+1/2 e−(z+r+1/2) 2
Z π/2 "√ #
r
2 v(θ) sin θ
× cos2z θ dθ (step II)
−π/2 log v(θ)

√
= (z + r + 1/2)z+1/2 e−(z+r+1/2) 2
Z π/2 " ∞
#
a0 (r) X
× cos2z θ + ak (r) cos (2kθ) dθ . (step III)
−π/2 2 k=0

The series in Equation (1.1) then follows from the last integral upon
integrating term by term.
The derivation is first retraced more or less along the same lines as
Lanczos uses in his paper [14]. This approach uses Fourier series tech-

40
Chapter 3. The Lanczos Formula

niques in a novel albeit complicated way to obtain explicit formulas

for the coefficients an (r) as a linear combination of Chebyshev polyno-
mial coefficients. The appearance of Chebyshev coefficients suggests a
connection with the orthogonal system of Chebyshev polynomials, and
indeed the derivation in the setting of these polynomials provides a
slightly cleaner argument. Finally, both of these derivations are seen to
be cases of a much more general idea, that of inner products in Hilbert
space.
The following notation which will be standard throughout the re-
mainder of this work:
Definition 3.1. Define
H0 (z) = 1 ,
and for k ≥ 1,
Γ(z + 1)Γ(z + 1)
Hk (z) =
Γ(z − k + 1)Γ(z + k + 1)

1
=
(z + 1)k (z + 1)−k

z · · · (z − k + 1)
= .
(z + 1) · · · (z + k)

Using this new notation, (1.1) may be restated1 :

∞
√ z+1/2 −(z+r+1/2)
X ′
Γ(z + 1) = 2π (z + r + 1/2) e ak (r)Hk (z) . (3.1)
k=0

3.1 The Lanczos Derivation

Lanczos uses as his starting point Euler’s formula (2.5) and transforms
the integral in a series of steps. The motivation for these transforma-
tions is not immediately obvious, though the effect is to extract from
1 The prime notation in the sum indicates that the k = 0 term is given weight
1/2.

41
Chapter 3. The Lanczos Formula

the integral a term similar in form to the first factor of Stirling’s series
(2.12). The integral which remains is then precisely the form required
to perform a clever series expansion.
In the exposition the following conventions have been adopted:

1. The complex variable z = σ + it where σ = Re(z) and t = Im(z);

2. For a, b ∈ C with a not a negative real number or zero, define

ab = eb log a where Im(log a) is between −π and π.

3.1.1 Preliminary Transformations

Beginning with
Z ∞
Γ(z + 1) = tz e−t dt, Re(z) > −1 , (3.2)
0

make the substitution t → αt where Re(α) > 0 to obtain

Z ∞
z+1
Γ(z + 1) = α tz e−αt dt .
0

The validity of this substitution for complex α can be shown using

contour arguments or via the following lemma2 :

Lemma 3.1. For Re(α) > 0 and Re(z) > −1,

Z ∞
z+1
Γ(z + 1) = α tz e−αt dt .
0

Proof of Lemma 3.1: For fixed z with Re(z) > −1 and α > 0 real,
replacing t with αt in (3.2) gives
Z ∞
z+1
Γ(z + 1) = α tz e−αt dt .
0

Viewed as a function of α, the right hand side is analytic in the right

half α-plane, and equals the constant Γ(z + 1) along the positive real
2 This result can be found as an exercise in [30, p.243]. For real α > 0 the
result is elementary, however Lanczos makes no comment as to its validity
for complex α.

42
Chapter 3. The Lanczos Formula

axis. Hence the integral must equal Γ(z + 1) throughout the right half
α-plane.

Next set α = z + r + 1, where r ≥ 0, so that

Z ∞
z+1
Γ(z + 1) = (z + r + 1) tz e−(z+r+1)t dt, Re(z) > −1 . (3.3)
0

Here r ≥ 0 is a free parameter which plays a key role later. Now reduce
the bounds of integration to a finite case via the change of variable
t = 1 − log v, dt = (−1/v)dv:
Z 0
z+1 −1
Γ(z + 1) = (z + r + 1) (1 − log v)z e−(z+r+1)(1−log v) dv
e v
Z e
z+1 −(z+r+1) 1
= (z + r + 1) e (1 − log v)z v z v r v dv
v
Z0 e
= (z + r + 1)z+1 e−(z+r+1) [v(1 − log v)]z v r dv . (3.4)
0

The integrand in this last expression is O(v σ+r−ǫ ) for every ǫ > 0 as
v → 0, while it is O(|v − e|σ ) as v → e . Equation (3.4) is therefore
valid for Re(z) > −1.3
The substitution t = 1 − log v appears to be the key to the Lanczos
method in the sense that it puts all the right pieces in all the right
places for the series expansion to come. It has the effect of peeling off
3 Thereis an error in [17, p.30, Eq.(1)] at this point concerning the domain
of convergence of the formula. There the author makes the replacement
z → z − 1/2 and states (using σ to signify what is here denoted r):

Γ(z + 1/2) = (z + σ + 1/2)z+1/2 exp [−(z + σ + 1/2)]F (z),

Re
F (z) = 0 [v(1 − log v)]z−1/2 v σ dv, Re(z + σ + 1/2) > 0 .

However, the integral F (z) diverges with, for example, z = −3/4 and a
choice of σ = 5/4. To see this, note that v(1 − log v) ≤ (e − v) on [0, e], so
that [v/(e − v)]5/4 ≤ [v(1 − log v)]−5/4 v 5/4 , whence
Z e Z e
5/4
[v/(e − v)] dv ≤ [v(1 − log v)]−5/4 v 5/4 dv .
0 0
Re
But the left hand side 0 [v/(e − v)]5/4 dv = ∞.

43
Chapter 3. The Lanczos Formula

leading terms similar in form to Stirling’s formula, while factoring the

integrand into separate powers of z and r. The sequence of substitutions
used here, however, namely α → (z + r + 1) and t → 1 − log v, is
quite different from the seemingly disconnected chain of substitutions
Lanczos uses in [14] to arrive at equation (3.4). His steps are perhaps
clues to the motivation behind his method and so these are reproduced
in Appendix A.

3.1.2 The Implicit Function v(θ)

The next step in the derivation is the transformation of the integral
in (3.4) into a form more amenable to Fourier methods. The path of
v(1 − log v) is shown in Figure 3.1 so that v may be implicitly defined
as a function of θ via

cos2 θ = v(1 − log v) ,

where θ = −π/2 corresponds to v = 0, θ = 0 to v = 1, and θ = π/2 to

v = e.

1
v(1 − log v)

0 1 e
v

Figure 3.1: Path of v(1 − log v)

Making the substitution v(1−log v) = cos2 θ, dv = 2 sin θ cos θ/ log v dθ

44
Chapter 3. The Lanczos Formula

in (3.4) then gives

π/2
2 sin θ cos θ
Z
z+1 −(z+r+1)
Γ(z + 1) = (z + r + 1) e (cos2 θ)z v r dθ
−π/2 log v
π/2
2v r sin θ
Z
= (z + r + 1)z+1 e−(z+r+1) cos2z+1 θ dθ . (3.5)
−π/2 log v

The integrand in (3.5) is O(|θ + π/2|2σ+2r+1 ) near θ = −π/2, O(1)

near θ = 0, but is O(|θ − π/2|2σ+1 ) as θ → π/2 so this expression in
once again valid for Re(z) > −1.
√
Moving a 2 outside the integral and replacing z with z − 1/2, the
reason for which will become apparent later4 , yields
Z π/2 "√ #
r
2v sin θ
Γ(z + 1/2) = Pr (z) cos2z θ dθ , (3.6)
−π/2 log v

where Pr (z) is defined

√
Pr (z) = 2 (z + r + 1/2)z+1/2 e−(z+r+1/2) .

Now denote by fr (θ) the term in square brackets in (3.6) and by fE,r (θ)
its even part [fr (θ) + fr (−θ)] /2. Noting that the odd part of the inte-
grand integrates out to zero finally puts (3.6) in the form
Z π/2
Γ(z + 1/2) = Pr (z) cos2z θ fE,r (θ) dθ (3.7)
−π/2

for Re(z) > −1/2 and r ≥ 0. Equation (3.7) is the starting point for
the series expansion in (3.1).
4 Shifting the argument to z −1/2 at the outset, that is, starting the derivation
with Z ∞
Γ(z + 1/2) = tz−1/2 e−t dt, Re(z) > −1/2
0
in place of equation (3.2), would seem more logical, considering the reason
for the shift is yet to come anyway. For some reason Lanczos does not do
this.

45
Chapter 3. The Lanczos Formula

3.1.3 Series Expansion of fE,r

The next step in the development is the replacement of fE,r (θ) by a
series which can be integrated term by term. Lanczos first constructs
a Taylor series in sin θ for fr (θ) but notes that the resulting series after
integrating is of slow convergence. Instead, he turns from the ‘extrapo-
lating’ Taylor series to an ‘interpolating’ series of orthogonal functions:
a Fourier series. The properties of fE,r will be examined later, but
for the time being assume that this function may be represented by a
uniformly convergent Fourier series on [−π/2, π/2]. Further, this series
will contain only cosine terms since fE,r is even. Thus
∞
1 X kπθ
fE,r (θ) = a0 (r) + ak (r) cos
2 π/2
k=1

∞
X ′
= ak (r) cos (2kθ) (3.8)
k=0

where the coefficients ak (r) are given by

π/2
2
Z
ak (r) = fE,r (θ) cos (2kθ) dθ .
π −π/2

On the surface this substitution would appear to complicate mat-

ters, or at least offer little improvement. It turns out, however, that
thanks to a handy identity, the resulting integrals evaluate explicitly to
the rational functions Hk (z) which exhibit the poles of Γ(z + 1). This
identity is stated as a lemma, the proof of which can be found in [17,
p.16]:

Lemma 3.2. For Re(z) > −1/2

Z π/2 √ Γ(z + 1/2)
cos2z θ cos (2kθ) dθ = π Hk (z) (3.9)
−π/2 Γ(z + 1)

Picking up where we left off with (3.7), and using (3.8) to replace

46
Chapter 3. The Lanczos Formula

fE,r (θ), we obtain

Z π/2 ∞
X ′
2z
Γ(z + 1/2) = Pr (z) cos θ ak (r) cos (2kθ) dθ
−π/2 k=0
∞ Z π/2
X′
= Pr (z) ak (r) cos2z θ cos (2kθ) dθ
k=0 −π/2

∞
√ Γ(z + 1/2) X′
= Pr (z) π ak (r)Hk (z) . (3.10)
Γ(z + 1) k=0

Now it is just a matter of eliminating Γ(z + 1/2) on both sides of

(3.10), the reason for replacing z with z − 1/2 in (3.5), and moving the
Γ(z + 1) to the left hand side to get the final form:
∞
z+1/2 −(z+r+1/2)
√ X ′
Γ(z + 1) = (z + r + 1/2) e 2π ak (r)Hk (z) , (3.11)
k=0

where again, Re(z) > −1/2, and r ≥ 0. We will see later that conver-
gence of this series can be extended to the left of Re(z) = −1/2.

3.2 Derivation in the Chebyshev Setting

The argument used to obtain (3.11) is now repeated, but this time
using the set of Chebyshev polynomials {Tk (x)} orthogonal on the in-
terval [−1, 1]. For a brief overview of Chebyshev polynomials and their
connection to Fourier expansions the reader is referred to Appendix B.
Beginning with equation (3.4),
Z e
Γ(z + 1) = (z + r + 1) z+1 −(z+r+1)
e [v(1 − log v)]z v r dv ,
0

this time define v(x) implicitly with

1 − x2 = v(1 − log v) , (3.12)
where x = −1 corresponds to v = 0, x = 0 to v = 1, and x = 1 to
v = e. This yields
Z 1
2xv r

z+1 −(z+r+1) 2 z
Γ(z + 1) = (z + r + 1) e (1 − x ) dx . (3.13)
−1 log v

47
Chapter 3. The Lanczos Formula

√
Once again let Pr (z) = 2 (z + r + 1/2)z+1/2 exp [−(z + √ r + 1/2)], re-
place z with z−1/2, and denote by fE,r (x) the even part of 2xv r / log v,
giving
Z 1
dx
Γ(z + 1/2) = Pr (z) (1 − x2 )z fE,r (x) √ . (3.14)
−1 1 − x2

The motivation this time for the replacement z → √ z−1/2 is a bit clearer:
with the introduction of the weight function 1/ 1 − x√2 , the integral
in (3.14) reveals itself as an inner product on L2[−1,1] (dx/ 1 − x2 ) which
has the Chebyshev polynomials as orthogonal basis.
Expanding fE,r (x) in a series of even Chebyshev polynomials (the
justification for which will follow) produces
∞
X ′
fE,r (x) = ck (r)T2k (x) (3.15)
k=0

so that
1 ∞
dx
Z X ′
2 z
Γ(z + 1/2) = Pr (z) (1 − x ) ck (r)T2k (x) √ ,
−1 k=0
1 − x2

which upon integrating term by term gives

∞
√ Γ(z + 1/2) X′
Γ(z + 1/2) = Pr (z) π ck (r)(−1)k Hk (z) .
Γ(z + 1) k=0

The term by term integration uses an equivalent form of the iden-

tity (3.9):
1 √ Γ(z + 1/2)
dx
Z
(1 − x2 )z T2k (x) √ = (−1)k π Hk (z) .
−1 1 − x2 Γ(z + 1)

Finally, multiplying through by Γ(z + 1)/Γ(z + 1/2) and writing

ak (r) = (−1)k ck (r) once again yields
∞
z+1/2 −(z+r+1/2)
√ X ′
Γ(z + 1) = (z + r + 1/2) e 2π ak (r)Hk (z) . (3.16)
k=0

48
Chapter 3. The Lanczos Formula

3.3 Derivation in the Hilbert Space Set-

ting
It is worth noting that both (3.7) and (3.14) show that Γ(z +1/2)/Pr (z)
manifests itself as an inner product in Hilbert space, so that (3.16) can
be deduced directly using Hilbert space theory. In fact, the requirement
that fE,r be equal to a uniformly convergent Fourier (resp. Chebyshev)
series is unnecessary, if we ask only that fE,r be square summable with
respect to Lebesgue measure.
Take for example the Fourier setting. Parseval’s theorem [24, p.91]
says that for even f, g ∈ L2 [−π/2, π/2] with g real,
π/2 ∞
2
Z
′
fˆk ĝk
X
f (θ)g(θ) dθ =
π −π/2 k=0

where fˆk denotes the Fourier coefficient

π/2
2
Z
fˆk = f (θ) cos (2kθ) dθ .
π −π/2

In the present case, the Fourier coefficients of p(θ) = cos2z (θ) are
π/2
2
Z
p̂k = cos2z (θ) cos (2kθ) dθ
π −π/2

2 Γ(z + 1/2)
=√ Hk (z) ,
π Γ(z + 1)

from Lemma 3.2, while those of fE,r are ak (r), assuming fE,r ∈ L2 [−π/2, π/2].

49
Chapter 3. The Lanczos Formula

Thus by Parseval’s theorem, (3.7) becomes simply

Z π/2
Γ(z + 1/2) = Pr (z) cos2z θ fE,r (θ) dθ
−π/2

∞
π X′
= Pr (z) p̂k ak (r)
2
k=0

∞
√ Γ(z + 1/2) X′
= Pr (z) π ak (r) Hk (z) ,
Γ(z + 1)
k=0

which again yields equation (3.1).

3.4 Additional Remarks

To conclude this chapter, a number of remarks on various aspects of
the derivation of (3.1) are worthy of mention.
The first remark is concerning the use of both Fourier and Cheby-
shev series to achieve the same result. These are clearly equivalent
since, as Boyd notes in [3], a Chebyshev polynomial expansion is merely
a Fourier cosine series in disguise. In this case, the transformation
x → sin θ links one to the other. Based on the extensive use of Cheby-
shev polynomials in his work [11], Lanczos was likely aware of the pos-
sibility of using either expansion in his derivation. Indeed, he does
make the substitution (3.12) in his paper, but not as a springboard
for a Chebyshev expansion, but rather to set the stage for expressing
f0 (θ) as a Taylor polynomial in sin θ. From the periodicity of this last
function he then argues that the interpolating Fourier cosine series is a
better choice in terms of convergence versus the extrapolating Taylor
series.
Secondly, the Fourier methods used here, in particular the use of
the identity (3.9), bear a striking resemblance to a technique used in
Lanczos’ 1961 work [13, p.45]. There he derives an expansion for Bessel
functions of arbitrary order in terms of Bessel functions of even integer
order. Part of that derivation involves a clever method for computing

50
Chapter 3. The Lanczos Formula

the coefficients of the expansion which will be described in the sequel.

This common thread is not so surprising since, as Lanczos notes in the
acknowledgments of each of [14] and [13], the winter of 1959-60 saw
much of both works completed at the Mathematics Research Center of
the U.S. Army at the University of Wisconsin in Madison.
Thirdly, the derivation of (3.1) can be simplified if one considers
first only real z > −1/2 followed by an application of the principle of
analytic continuation at the end. This approach obviates Lemma 3.1
but requires consideration of the analytic properties of the infinite series
factor of the formula. Such matters will be examined later when (3.1)
is extended to the left of Re(z) = −1/2.
The final remark involves the free parameter r as it relates to the
Fourier coefficients ak (r) (or equivalently ck (r) in the Chebyshev set-
ting). There are several methods for computing the coefficients and
these will be examined in due course. For the moment, note Lanczos’
own formulation: the constant term of the series is
1/2
2e
a0 (r) = er , (3.17)
π(r + 1/2)

while for k ≥ 1,
k
2X
ak (r) = C2j,2k Fr (j) . (3.18)
π j=0

Here Fr (j) = 2−1/2 Γ(j + 1/2)(j + r + 1/2)−j−1/2 exp (j + r + 1/2) and

C2j,2k is the coefficient of x2j in the 2k th Chebyshev polynomial.
These coefficients are functions of r and this parameter plays a
fundamental role in the rate of convergence of the series. This role
will be examined in detail, but for now observe that the constraint
r ≥ 0 imposed by Lanczos is unduly restrictive in the sense that the
substitution in equation (3.3) requires only that Re(r) ≥ 0. Under this
relaxed constraint, the Fourier coefficients given by (3.17) and (3.18) are
then single valued analytic functions of r in this region. The properties
of the ak (r) as a function of a complex variable were not considered
here, but it would be interesting to investigate further the effect of this
generalization on (3.1).

51
Chapter 4
The Functions v and fE,r

The derivation of Lanczos’ main formula (1.1) in Chapter 3 relies on

the validity of expanding the implicitly defined function fE,r (θ) of equa-
tion (3.7) (resp. fE,r (x) of (3.14)) into an infinite interpolating Fourier
(resp. Chebyshev) series. Furthermore, the smoothness and summa-
bility properties of fE,r determine the asymptotic growth rate of the
coefficients ak (r) (resp. ck (r)), and consequently the convergence of
the series (1.2).
In this chapter the properties of the implicit functions v, fr and fE,r
are examined in detail, and bounds are established on the coefficients
appearing in the series (3.8) (resp. (3.15)). The principal results are

(i) In the Chebyshev setting, fE,r (x) ∈ C ⌊r⌋ [−1, 1] and is of bounded
variation, so that the expansion (3.15) is justified.

(ii) In the Fourier setting, fE,r (θ) ∈ C ⌊2r⌋ [−π/2, π/2] and is of bounded
variation, thus justifying the expansion (3.8).
(n)
(iii) fE,r ∈ L1 [−π/2, π/2], where n = ⌈2r⌉.

(iv) ak (r) = O(k −⌈2r⌉ ) as k → ∞.

4.1 Closed Form Formulas

We begin by finding explicit expressions for v and fE,r in terms of
Lambert W functions W0 and W−1 . This is useful for several reasons:
the first is that graphs of v and fE,r can be easily plotted using the built
in Lambert W evaluation routines of Maple 8. The second is that the

52
Chapter 4. The Functions v and fE,r

smoothness properties of v and fE,r can be deduced directly from those

of W0 and W−1 . Third, the expression for fE,r can be used to compute
approximate values of the coefficients ak (r) using finite Fourier series.
In this section we work first in the Chebyshev setting and establish
explicit formulas for v(x) and fE,r (x) on [−1, 1]. These formulas are in
terms of the two real branches W−1 and W0 of the Lambert W function
defined implicitly by
W (t)eW (t) = t .
See [5] for a thorough treatment of Lambert W functions.

4.1.1 Closed Form for v(x)

Recall that v(x) was defined implicitly via

1 − x2 = v(1 − log v) , (4.1)

where x = −1 corresponds to v = 0, x = 0 to v = 1, and x = 1 to

v = e. Letting w = log (v/e), this is equivalent to

x2 − 1
= wew
e
whence
x2 − 1

w=W
e
so that
2
v x −1
= exp W
e e
2
x −1
e
= x2 −1
. (4.2)
W e

For −1 < x < 0, 0 < v < 1 which corresponds to the real branch
W−1 of the Lambert W function. For 0 < x < 1, 0 < v < e which
corresponds to the branch W0 . Using the Heaviside function H(x), and
letting y(x) = (x2 − 1)/e, v can thus be expressed as the single formula

v(x) y(x) y(x)

= H(−x) + H(x) . (4.3)
e W−1 (y(x)) W0 (y(x))

53
Chapter 4. The Functions v and fE,r

In developing this expression, care must be taken to ensure that v takes

on the prescribed values at the endpoints x = 1, −1 where W−1 (y(x)) →
−∞, and also at x = 0 where H(x) has a jump discontinuity. Taking
limits of (4.3) at these points shows that indeed v takes on the pre-
scribed values and is thus continuous on [−1, 1]. See Figure 4.1 for a
plot of v(x)/e.

v(x)
e

−1 0 1

Figure 4.1: v(x)/e, −1 ≤ x ≤ 1

4.1.2 Closed Form for fr (x) and fE,r (x)

√
For the function fr (x) = 2 xv r / log v, note that log v = 1 + W ((x2 −
1)/e) from (4.2), so that fr (x) may be expressed
 y(x) r
y(x)
r 
fr (x) √ W−1 (y(x)) W0 (y(x))
r
= 2x H(−x) + H(x) ,
e 1 + W−1 (y(x)) 1 + W0 (y(x))

54
Chapter 4. The Functions v and fE,r

while its even part fE,r (x) = [fr (x) + fr (−x)] /2 becomes
 y(x) r
y(x)
r 
fE,r (x) |x|  W0 (y(x)) W−1 (y(x))
= √ −  . (4.4)
er 2 1 + W0 (y(x)) 1 + W−1 (y(x))

√
Plots of 2 fE,r (x)/er are shown in Figure 4.2 for various values of
r. Notice the value of the function at the endpoints is one√for the four
values√of r plotted. It is clear from the definition fr (x) = 2 xv r / log v
that 2 fE,r (±1)/er = 1 if r ≥ 0 and is +∞ otherwise. Also notice
the extreme behaviour of the derivative of fE,r at the end points for
r = 0. Lanczos notes that the parameter r serves to smooth out this
singularity and hence improve the convergence of the Fourier series of
fE,r .

r=0

r=1

r=2

r=4

−1 0 1
√
Figure 4.2: 2 fE,r (x)/er , −1 ≤ x ≤ 1, r = 0, 1, 2, 4

55
Chapter 4. The Functions v and fE,r

4.2 Smoothness of v r (x), fr (x) and fE,r (x)

The smoothness of v r , fr and fE,r can √be seen to be essentially the
same problem with the observation that 2 fr (x) = d/dx[v r+1/(r + 1)].
This follows from the fact that 1 − x2 = v(1 − log v) from which

d 2x
v(x) = . (4.5)
dx log v

In this section it is shown that v r , fr (and hence fE,r ) have ⌊r⌋ contin-
uous derivatives on [−1, 1].

Smoothness on (−1, 1]

Away from x = −1 (where v > 0), the smoothness of v r is precisely the

same as that of v. As such, it is enough to consider v(x) on (−1, 1].
On (−1, 0) ∪ (0, 1), the functions Wk ((x2 − 1)/e) of (4.3) are smooth
and not zero (see [5]), and so v(x) is smooth on this set. About x =
1, v(x)/e = exp W0 ((x2 − 1)/e) which is again smooth since W0 (t) is
smooth at t = 0.
This leaves x = 0 where the situation is not so clear, but the implicit
function theorem gives the result. By expressing the right hand side
of (4.1) as a Taylor series, it is not difficult to show that (4.1) may be
written "∞ #1/2
X
j (v − 1)j
x = (v − 1) (−1)
j=0
(j + 1)(j + 2)

where the series on the right hand side converges absolutely and uni-
formly for |v − 1| < 1. Denote this right hand side √ by h(v) and set
φ(v, x) = x−h(v). Now φ(1, 0) = 0, ∂φ/∂v|(1,0) = −1/ 2, and φ ∈ C ∞
at (1, 0). Therefore by the implicit function theorem there are open sets
U ⊂ R2 and W ⊂ R and a function g : W → R such that g ∈ C ∞ ,
(1, 0) ∈ U, 0 ∈ W , g(0) = 1, and φ(g(x), x) = 0 for every x ∈ W (see
[19, p.103]). Thus v = g(x) is smooth about x = 0.
Thus v is smooth as a function of x on (−1, 1], from which v r and
fr are as well.

56
Chapter 4. The Functions v and fE,r

Smoothness at x = −1

It remains to examine smoothness at x = −1. At this point the role

of the parameter r becomes crucial and must be taken into account.
When r = 1, v and v ′ = 2x/ log v are zero as x → −1+ , while v ′′ =
2/ log v − 4x2 /[v(log v)3 ] → ∞ there. Thus v is once continuously
differentiable at x = −1. Raising v to the power r > 1, on the other
hand, forces v r to zero more rapidly as x =→ −1+ , and so has the
effect of smoothing out the singularity there.
Begin with a simple calculus result: if r ≥ 0 and k ≤ 0 then

lim + xj v r (log v)k = lim+ v r |log v|k

x→−1 v→0
< ∞. (4.6)

With this in mind, consider the smoothness properties of the gen-

eralized function
g(x) = xj v r (log v)k
where j and k are integers and r is real. Using (4.5) to differentiate
g(x) yields

d j r
x v (log v)k =

dx
jxj−1 v r (log v)k + 2rxj+1 v r−1 (log v)k−1 + 2kxj+1 v r−1 (log v)k−2
which can be abstracted as
d/dx
G(j, r, k) −→ jG(j−1, r, k)+2rG(j+1, r−1, k−1)+2kG(j+1, r−1, k−2) .
(4.7)
In the present case, v r corresponds to G(0, r, 0), while fr corresponds
to 2G(1, r, −1). Now beginning with G(0, r, 0) and differentiating re-
peatedly, (4.7) shows that with each derivative, the k component of
any resulting G terms remains zero or negative while the least r com-
ponent reduces by 1. By (4.6), at least the first ⌊r⌋ derivatives of v r
are bounded as x =→ −1+ . The case r = 1 shows that this bound is
best possible. Since fr (x) = d/dx[v r+1 /(r + 1)], it follows that fr has
⌊r⌋ derivatives as well.
In summary

57
Chapter 4. The Functions v and fE,r

Theorem 4.1.

(i) fr (x) and v r are smooth on (−1, 1];

(ii) fr (x), fE,r (x), and v r are in C ⌊r⌋ [−1, 1].

Before proceeding to the examination of fr as a function of θ, we

end this section with a result which justifies the representation of fr in
equation (3.15) by a uniformly convergent Chebyshev series:

Theorem 4.2. fr (x) is increasing on [−1, 1].

Proof of Theorem 4.2: It is enough to show fr′ (x) ≥ 0 on (−1, 1).

Since √
2 xv r
fr (x) =
log v
and v r is non negative and increasing, it is enough to show that
x
h(x) =
log v

is non negative and increasing. It is clear that h(x) = v ′ (x)/2 ≥ 0.

Now
log v − 2x2 /(v log v)
h′ (x) = ,
(log v)2
so it is enough to show that the numerator

p(x) = log v − 2x2 /(v log v)

v (log (v))2 − 2 + 2 v − 2 v log (v)
=
v log (v)

is non negative. But this last function is minimized when v = 1 where

p(0) = 0, completing the proof.

Since fr (x) is increasing its total variation is fr (1) − fr (−1), from

which

58
Chapter 4. The Functions v and fE,r

Corollary 4.1. fr (x) and fE,r (x) are functions of bounded variation
on [−1, 1].

From the continuity of fE,r (x) and Corollary 4.1 it follows that
fE,r (x) may be expressed as a uniformly convergent Chebyshev series
on [−1, 1].

4.3 Smoothness of v r (θ), fr (θ) and fE,r (θ)

By setting x = sin θ and
d dx d
= ,
dθ dθ dx
it follows easily from the previous section that v r (θ), fr (θ) and fE,r (θ)
have at least ⌊r⌋ continuous derivatives on [−π/2, π/2]. In fact, this
argument shows that these functions are smooth on (−π/2, π/2], with
a singularity at θ = −π/2.
In this section, the precise order of the singularity of fr (θ) (and
consequently fE,r (θ)) at θ = −π/2 is determined. This will then be
used to bound the growth of the Fourier coefficients ak (r). To do so, we
proceed as in the Chebyshev setting and define a generalized function,
although the generalization is slightly more complicated in this case.

4.3.1 Generalized fr (θ)

Define G(β, j, k, l; θ) as

G(β, j, k, l; θ) = v β sinj θ cosk θ(log v)l , (4.8)

√
where β ∈ R, j, k, l ∈ Z. Note that fr (θ) = 2 G(r, 1, 0, −1; θ). From
tedious application of derivative rules and using the fact that
2 sin θ cos θ
v ′ (θ) = ,
log v

59
Chapter 4. The Functions v and fE,r

G(β, j, k, l; θ) has the derivative form

d
[G(β, j, k, l; θ)] = 2βG(β − 1, j + 1, k + 1, l − 1; θ)
dθ
+jG(β, j − 1, k + 1, l; θ)
−kG(β, j + 1, k − 1, l; θ)
+2lG(β − 1, j + 1, k + 1, l − 2; θ) . (4.9)

Thus G and its derivatives can always be written as a linear combination

of functions of the form (4.8). The following lemmas will help to clarify
the dependence of G on β, j, k and l. Begin with an easy result from
calculus:

Lemma 4.1. Suppose β ∈ R and l ∈ Z. Then



 0 if β >0
 ∞ if β <0


β l
lim v |log v| = 1 if β = 0, l = 0
v→0+
∞ if β = 0, l > 0




0 if β = 0, l < 0 .


From this follows

Lemma 4.2. Suppose β ∈ R and j, k, l ∈ Z. Then



 0 if β + k/2 > 0
l+j
 (−1) · ∞ if β + k/2 < 0


lim + G(β, j, k, l; θ) = (−1)l+j if β + k/2 = 0 and k/2 + l = 0
π
θ→− 2 l+j
(−1) · ∞ if β + k/2 = 0 and k/2 + l > 0




0 if β + k/2 = 0 and k/2 + l < 0 .


Proof of Lemma 4.2: Near θ = −π/2

cos θ = v 1/2 |1 − log v|1/2 .

Also recall that as θ → −π/2+ , v → 0+ . Now write

60
Chapter 4. The Functions v and fE,r

lim + G(β, j, k, l; θ) = lim + v β sinj θ cosk θ(log v)l

θ→− π
2
θ→− π
2

= lim v β v k/2 |1 − log v|k/2 sinj θ |log v|l (−1)l

+
θ→− π
2

= lim + v β+k/2 |log v|l+k/2 (−1)l+j

θ→− π
2

and use Lemma (4.1).

The next section shows how to use these lemmas to determine lower
bounds on the order to which G is continuously differentiable at −π/2.

4.3.2 Derivatives of G(β, j, k, l; θ)

To apply the results of the previous section to fr , abstract the process
of differentiating G(β, j, k, l; θ) as follows. Identify with G(β, j, k, l; θ)
the monomial xβ y j z k w l , and view differentiation of G(β, j, k, l; θ) as the
linear mapping :

D : xβ y j z k w l 7→ 2βxβ−1 y j+1z k+1 w l−1

+jxβ y j−1z k+1 w l
−kxβ y j+1z k−1 w l
+2lxβ−1 y j+1z k+1 w l−2 . (4.10)

Higher order derivatives of G(β, j, k, l; θ) are then obtained recur-

sively as D2 = D ◦ D, D 3 = D ◦ D ◦ D, etc. For consistency, let D 0
denote the identity mapping.

4.3.3 Application to fr
We arrive finally at the analysis of fr .

Theorem 4.3. fr (θ) is continuous, increasing and hence of bounded

variation on [−π/2, π/2].

61
Chapter 4. The Functions v and fE,r

Proof of Theorem 4.3: This follows from Theorems 4.2 and 4.1 upon
writing x = sin θ and
dfr dfr dx
= .
dθ dx dθ

Theorem 4.4. fr (θ) has ⌊2r⌋ continuous derivatives at θ = −π/2.

Proof√of Theorem 4.4: Using the notation of Section (4.3.1), write

fr (θ)/ 2 = G(r, 1, 0, −1; θ). By Lemma 4.2, fr is continuous at θ =
−π/2 for r ≥ 0. Now consider the associated sequence

D0 (xr y 1 z 0 w −1 ), D1 (xr y 1 z 0 w −1 ), D2 (xr y 1z 0 w −1 ), . . . (4.11)

The task is to determine the least n such that D n+1 (xr y 1z 0 w −1 )

contains a term whose powers meet one of the ∞ conditions in Lemma
4.2. G(β, j, k, l; θ) will then be n-times continuously differentiable at
−π/2.
Observe first that the third and fourth conditions of Lemma 4.2
(β + k/2 = 0 and k/2 + l ≥ 0) can never be satisfied by a term of one of
the D n (xr y 1z 0 w −1 ) in the sequence (4.11), since these conditions imply
β − l ≤ 0, which translates into deg (x) − deg (w) ≤ 0 in some term of
D n (xr y 1 z 0 w −1 ). But deg (x) − deg (w) = r + 1 > 0 in D 0 and (4.10)
implies this quantity is non-decreasing with n.
The only other ∞ condition in Lemma 4.2 is satisfied when deg (x)+
deg (z)/2 < 0 for a term in some D n (xr y 1z 0 w −1 ). In D 0 , deg (x) +
deg (z)/2 = r ≥ 0. Equation (4.10) implies that the least value
of deg (x) + deg (z)/2 taken over the summands of D n (xr y 1z 0 w −1 )
decreases by at most 1/2 with each increase in n. Further, since
deg (w) = −1 in D0 , the fourth term of (4.10) guarantees that this
decrease will indeed occur with each n. Thus, it is enough to determine
the largest integer n such that r − n(1/2) ≥ 0, which is the largest n
such that n ≤ 2r, i.e. n = ⌊2r⌋. In other words, for r ≥ 0, fr (θ) has
⌊2r⌋ continuous derivatives. Lemma 4.2 also gives that all of these ⌊2r⌋
derivatives are zero.

(n)
Theorem 4.5. fr (π/2) = 0 for n a positive odd integer.

62
Chapter 4. The Functions v and fE,r

√
Proof of Theorem 4.5: Again identify fr (θ)/ 2 with G(r, 1, 0, −1; θ)
and consider the associated sequence of derivatives

D 0 (xr y 1z 0 w −1), D1 (xr y 1 z 0 w −1 ), D2 (xr y 1 z 0 w −1 ), . . .

By (4.10), each term of every odd order derivative has deg (z) > 0. Thus
each term of every odd order derivative has a factor of cos θ which is
zero at θ = π/2.

(n)
Theorem 4.6. fr ∈ L1 [−π/2, π/2] where n = ⌈2r⌉.

Proof of Theorem 4.6: It has been established that f is smooth on

(n)
(−π/2, π/2], so it is sufficient to bound the growth offr as θ → −π/2+ .
Again it is easiest to work with the generalized fr as given in Section
4.3.1, from which it is not difficult to determine conditions on β, j, k
and l which guarantee summability of the integral. Beginning with
Z π/2
I = |G(β, j, k, l; θ)| dθ
−π/2
Z π/2 β j
v sin θ cosk θ

= dθ ,
−π/2
(log v)l

change to v variables with the substitution v(1 − log v) = cos2 θ:

Z e β
v [1 − v(1 − log v)]j/2 [v(1 − log v)]k/2 log v

I = (log v)l [1 − v(1 − log v)]1/2 [v(1 − log v)]1/2 dv

0
Z e
= v β [1 − v(1 − log v)](j−1)/2 [v(1 − log v)](k−1)/2 |log v|1−l dv
0

As v → 0+ , the integrand is O(v β+k/2−1/2 (log v)k/2−l+1/2 ). Provided

β + k/2 − 1/2 > −1, that is, β + k/2 + 1/2 > 0, the integral converges.
√
In the case at hand of fr (θ) = 2 G(r, 1, 0, −1; θ), the exponent
β + k/2 − 1/2 starts off at r − 1/2 and decreases by 1/2 with each
derivative as (4.10) shows. Thus if n is chosen to be the largest integer
such that
r − 1/2 − n(1/2) > −1

63
Chapter 4. The Functions v and fE,r

R π/2 (n)
then convergence of −π/2
|fr | dθ is guaranteed. That is, if n is the
(n)
largest integer such that n < 2r + 1, which is just ⌈2r⌉, then fr ∈
L1 [−π/2, π/2].

The results of Theorems 4.3, 4.4, 4.5 and 4.6 extend easily to the
even part of fr :

Corollary 4.2. Suppose r ≥ 0 and

fr (θ) + fr (−θ)
fE,r (θ) = .
2
Then

(i) fE,r is a continuous function of bounded√variation on [−π/2, π/2]

such that fE,r (−π/2) = fE,r (π/2) = er / 2 and fE,r (0) = 1;
(n) (n)
(ii) fE,r ∈ C n [−π/2, π/2], where n = ⌊2r⌋, and fE,r (−π/2) = fE,r (π/2) =
0 for n ≤ ⌊2r⌋ an odd integer;
(n)
(iii) fE,r ∈ L1 [−π/2, π/2], where n = ⌈2r⌉.

4.4 Fourier Series and Growth Order of

Coefficients
We conclude this chapter with the application of Corollary 4.2 to the
justification of the replacement of fE,r by its Fourier series in equa-
tion (3.8), and to the estimate of the growth order of the Fourier coef-
ficients ak (r) in equation (1.1).
First of all, from (i) of Corollary 4.2 it follows that fE,r has a Fourier
series expansion which converges uniformly on [−π/2, π/2] (see [25]).
Since fE,r (θ) is even, this series will contain only cosine terms. Thus
∞
1 X
fE,r (θ) = a0 (r) + ak (r) cos (2kθ)
2 k=1

64
Chapter 4. The Functions v and fE,r

for coefficients
π/2
2
Z
ak (r) = fE,r (θ) cos (2kθ) dθ .
π −π/2

By substituting θ = π/2 and θ = 0 at this stage we get

∞
er 1 X
√ = a0 (r) + (−1)k ak (r)
2 2 k=1

and ∞
1 X
1 = a0 (r) + ak (r) , (4.12)
2 k=1

useful facts for testing the accuracy of computed ak (r)’s and examining
the error later.
Conclusions (ii) and (iii) of Corollary 4.2 allow the determination
of a growth bound on the ak (r):
Theorem 4.7. As k → ∞, ak (r) = O(k −⌈2r⌉ ).

Proof of Theorem 4.7: Beginning with

2 π/2
Z
ak (r) = fE,r (θ) cos (2kθ) dθ ,
π −π/2
integrate by parts n = ⌈2r⌉ times to get
(−1)⌈r⌉ 2 π/2 (n)
Z
ak (r) = f (θ) {cos (2kθ) | sin (2kθ)}n dθ . (4.13)
(2k)n π −π/2 E,r
Here

cos (2kθ) if ⌈2r⌉ is even,
{cos (2kθ) | sin (2kθ)}n =
sin (2kθ) if ⌈2r⌉ is odd.
From (4.13) it follows that
π/2
1 2
Z
(n)
|ak (r)| ≤ fE,r (θ) dθ ,
(2k)n π −π/2

whence
ak (r) = O(k −⌈2r⌉ ) . (4.14)

65
Chapter 5
Analytic Continuation of
Main Formula

At the end of the derivation of the main formula (1.1) in Section 3.1, it
X′
was noted that convergence of the series Sr (z) = ak (r)Hk (z) could
be extended to the left of the line Re(z) = −1/2. In this chapter this
statement is made precise by proving in detail that the series converges
absolutely and uniformly compact subsets of the half plane Re(z) >
−r excluding the negative integers. It is then shown that the main
formula (1.1) is analytic, and thus defines Γ(z + 1) in this region.
Begin with some notation:
Definition 5.1. For r ≥ 0, define the open set

Ωr = {z ∈ C | Re(z) > −r and z 6= −1, −2, −3, . . .} .

5.1 Domain of Convergence of Sr (z)

Theorem 5.1. Suppose r ≥ 0. Then the series
∞
X ′
Sr (z) = ak (r)Hk (z)
k=0

converges absolutely and uniformly on every compact K ⊂ Ωr . Conse-

quently, the series defines an analytic function on Ωr .

Proof of Theorem 5.1: First, observe that each term of the series is
analytic on Ωr . Now suppose K is a compact subset of Ωr , let D(0, R)

66
Chapter 5. Analytic Continuation of Main Formula

be a closed disk of radius R containing K, and let δ be the distance

between K and C \ Ωr .
It will first be shown that |Hk (z)| = O k 2r−2δ−1 on K. For k > R,

by the reflection formula (2.10),

Γ(z + 1)Γ(z + 1)
Hk (z) =
Γ(z + 1 + k)Γ(z + 1 − k)

(−1)k+1 [Γ(z + 1)]2 sin (πz) Γ(1 − z + k)

=
π Γ(1 + z + k)(k − z)

Γ(1 − z + k)
= h(z) , say.
Γ(1 + z + k)(k − z)

Now h(z) is analytic on Ωr , so |h(z)| is bounded on K since K is

compact. It remains to bound the second factor. For k > R, 1 − z + k
and 1 + z + k lie in the right half plane Re(z) ≥ 0, so using Stirling’s
formula to replace the gamma functions in the fraction,

Γ(1 − z + k)
Γ(1 + z + k)(k − z)

√
2πez−k (k − z)k−z+1/2 (1 + δ1 (k)) 1
=√
2πe−z−k (k + z)k+z+1/2 (1 + δ2 (k)) (k − z)

k z 1/2
2z k−z 1 k−z 1 (1 + δ1 (k))
=e
k+z (k − z)(k + z) k+z (k − z) (1 + δ2 (k))

where δj (k) → 0 as k → ∞ and the δj do not depend on z. Selecting

67
Chapter 5. Analytic Continuation of Main Formula

z ∈ D(0, R) to make this last expression as large as possible,

Γ(1 − z + k)

Γ(1 + z + k)(k − z)

k 1/2
2R k+R k+R (1 + δ1 (k)) 1
≤e
k−R k−R (1 + δ2 (k)) (k − R)2σ+1

k 1/2
2R k+R k+R (1 + δ1 (k)) 1
≤e
k−R k−R (1 + δ2 (k)) (k − R)2(−r+δ)+1

= O k 2r−2δ−1

.

This last inequality follows since, as k → ∞,

k 1/2
2R k+R k+R (1 + δ1 (k))
e → e4R .
k−R k−R (1 + δ2 (k))

Recalling that |h(z)| is bounded, the conclusion is

|Hk (z)| = O k 2r−2δ−1

on K.
Since |ak (r)| = O(k −⌈2r⌉ ) from Theorem 4.7, it then follows that

|ak (r)Hk (z)| = O(k 2r−⌈2r⌉−1−2δ )

= O(k −1−2δ ) ,

whence, for some positive constant A and any n ≥ 1

∞ ∞
X X 1
|ak (r)Hk (z)| ≤ A <∞,
k=n k=n
k 1+2δ

completing the proof.

68
Chapter 5. Analytic Continuation of Main Formula

5.2 Analytic Continuation

Theorem 5.2. For r ≥ 0, the function
∞
√ z+1/2 −(z+r+1/2)
X ′
2π (z + r + 1/2) e ak (r)Hk (z)
k=0

is analytic on Ωr .
√
Proof of Theorem 5.2: The factor 2π (z + r + 1/2)z+1/2 is analytic
since Re(z + r + 1/2) > 0 on Ωr . The factor exp [−(z + r + 1/2)] is
entire. The series factor is analytic by Theorem 5.1.

Theorem 5.3.
∞
√ z+1/2 −(z+r+1/2)
X ′
Γ(z + 1) = 2π (z + r + 1/2) e ak (r)Hk (z)
k=0

on Ωr .

Proof of Theorem 5.3: The right hand side is analytic on Ωr and

equals Γ(z + 1) on Re(z) > −1/2. Hence by the principle of analytic
continuation, it equals Γ(z + 1) throughout Ωr .

5.3 Additional Remarks

In his original paper, in reference to the main formula (1.1), Lanczos
remarks simply:

It is of interest to observe, however, that the convergence

of the infinite expansion extends even to the negative realm
and is in fact limited by the straight line

Re(z) > −(r + 1/2) .

69
Chapter 5. Analytic Continuation of Main Formula

This is most certainly true. Indeed, with r = 0 this bound is in agree-

ment with the bound given following equation (3.11). However, the
proof of Lanczos’ bound eludes me.
Observe that Sr (z) may be expressed

Sr (z) = Γ(z + 1)(2π)−1/2 (z + r + 1/2)−(z+1/2) ez+r+1/2 ,

and that, except for the simple poles at z = −1, −2, . . . which are
present on both sides of the equation, z = −r −1/2 is the first singular-
ity encountered as z moves into the left-hand plane. Lanczos’ statement
seems to be based on a principle similar to that for the power series
of an analytic function, that the radius of convergence of the series
extends to the first singularity. Closer in spirit is a result from the the-
ory of Dirichlet series: suppose f (s) is an analytic function P on the half
plane Re(s) > c1 , and f (s) is expressible as a series f (s) = ∞ k=1 qk k
−s

on Re(s) ≥ c2 ≥ c1 , where the coefficients qk are eventually positive.

Then the series actually converges on the larger set Re(s) > c1 .
The best that can be achieved with the techniques used here is
convergence on Re(z) > −⌈2r⌉/2, a slight improvement over Theo-
rem 5.3. If the bound ak (r) = O(k −⌈2r⌉ ) could be improved to ak (r) =
O(k −(2r+1) ) then the theorem would give absolute and uniform conver-
gence of Sr (z) on compact subsets of Re(z) > −(r + 1/2). Numerical
checks do support this slightly more rapid rate of decrease. For exam-
ple, in Figure 5.1, − log |ak (r)k 2r+1 | is compared against − log |ak (r)k ⌈2r⌉ |,
for the first 150 coefficients with r = 1.25 (− log |ak (r)| is also plotted
for reference). Of the two lower curves, the flatter one corresponding
to − log |ak (r)k 2r+1 | indicates the tighter bound.
The difficulty lies in finding a more precise bound on the Fourier
coefficients than the rather cumbersome estimate ak (r) = O(k −⌈2r⌉ ).
To do this requires a better understanding of the asymptotic behaviour
of fr (θ) about the singularity θ = −π/2. Letting h = (θ + π/2)2 , it is
true that
√ hr
fr (θ) ∼ 2
(− log h)r+1
as θ → −π/2. These two functions have the same number of derivatives
there, so the asymptotic tendency of their Fourier coefficients should
be the same. The asymptotic behaviour of the Fourier coefficients of
hr /(− log h)r+1 does not appear straightforward, however.

70
Chapter 5. Analytic Continuation of Main Formula

25

− log |ak (r)|

− log |ak (r)k ⌈2r⌉ |

− log |ak (r)k 2r+1 |

0
0 k 150
−(2r+1)
Figure 5.1: Growth of ak (r), k , and k −⌈2r⌉ with r = 1.25

71
Chapter 6
Computing the Coefficients

In order to use the formula (1.1) in practice, the Fourier coefficients

ak (r) must first be determined. The method given by Fourier theory,
namely the direct integration
2 π/2
Z
ak (r) = fE,r (θ) cos (2kθ) dθ (6.1)
π −π/2
in not practical due to the complicated nature of the factor fE,r (θ)
in the integrand. There are, however, a number of other methods for
computing the ak (r), and in this chapter these various methods are
discussed along with some comments on their practicality.
The first method is that of Lanczos as described in Chapter 1. The
second method is a recursive procedure which uses nested calculations
similar to Horner’s rule for increased efficiency. This method was found
to work very well in practice, and was used for most large scale calcula-
tions and numerical investigations. These first two methods rely on the
fact that Γ(z + 1) is explicitly known at the integers and half integers.
The third method is another recursive procedure based on the identity
Γ(z + 1) = zΓ(z). The fourth method uses the Chebyshev version of a
finite Fourier transform on the closed form expression for fE,r from Sec-
tion 4.1.2. This last method is essentially a numerical approximation
of the integral (6.1).

6.1 The Lanczos Method

The first method for finding the ak (r) is the clever and elegant method
of Lanczos, as given in his original paper, but carried out in the Cheby-
shev setting which simplifies it slightly. (Recall that ak (r) = (−1)k ck (r).)

72
Chapter 6. Computing the Coefficients

First define Fr (z) as the integral function

1
dx
Z
Fr (z) = (1 − x2 )z fE,r (x) √
−1 1 − x2

appearing in equation (3.14), and notice by this same equation that Fr

can also be written

Fr (z) = 2−1/2 Γ(z + 1/2)(z + r + 1/2)−z−1/2 exp (z + r + 1/2) .

Now denote the 2k th Chebyshev polynomial by T2k (x) = kj=0 C2j,2k x2j ,
P
and from (B.2) the coefficients in the Chebyshev series are defined as

2 1 dx
Z
ck (r) = fE,r (x)T2k (x) √ .
π −1 1 − x2
The Chebyshev series coefficients are then given by

2 1 dx
Z
ck (r) = fE,r (x)T2k (x) √
π −1 1 − x2
" k #
2 1 dx
Z X
= fE,r (x) C2j,2k x2j √
π −1 j=0
1 − x2
Z 1 " k
#
2 X dx
= fE,r (x) (−1)k C2j,2k (1 − x2 )j √ using (B.3)
π −1 j=0
1 − x2
k 1
2X dx
Z
k
= (−1) C2j,2k (1 − x2 )j fE,r (x) √
π j=0 −1 1 − x2
k
2Xk
= (−1) C2j,2k Fr (j) . (6.2)
π j=0

With a ready list of Chebyshev polynomial coefficients and precom-

puted Fr values, (6.2) can be concisely expressed
    
C0,0 Fr (0) c0 (r)
 C0,2 C2,2
2   Fr (1)   −c1 (r) 
 ..   ..  =  ..  .
   
π .  .   . 
k
C0,2k C2,2k · · · C2k,2k Fr (k) (−1) ck (r)

73
Chapter 6. Computing the Coefficients

In other words, ck (r) is simply a weighted sum of Fr values over

the first k integers, where the weights are Chebyshev polynomial co-
efficients. This is yet another peculiar feature of Lanczos’ paper: the
coefficients of the approximating series Sr (z) are in terms of the very
function being approximated. From this point of view, Sr (z) is the infi-
nite interpolation of (2π)−1/2 Γ(z+1)(z+r+1/2)−z−1/2 exp (z + r + 1/2)
by rational functions Hk (z).
To get an explicit form for the coefficients, an expression for the
coefficients C2j,2k is required. For j = 0, C0,2k = (−1)k , while for j ≥ 1

k−j k k+j
C2j,2k = (−1) 4j .
k+j k−j
p
Using this in equation (6.2) yields c0 (r) = 2e/[π(r + 1/2)] er and
r k j+1/2
2 r X j (k + j − 1)! e
ck (r) = ek (−1) .
π j=0
(k − j)!j! j + r + 1/2

The corresponding ak (r) appearing in (1.1) are then ak (r) = (−1)k ck (r).
Although this method provides a convenient closed form for com-
puting individual coefficients, it has the drawback of requiring generous
floating point precision to handle the addition of large values of alter-
nating sign which arise in intermediate calculations. This problem can
be partially overcome with the observation that each summand of (6.2)
contains a factor of exp (r), so this term can be factored out to reduce
the size of intermediate calculations. Even so, the problem persists.
For example, consider the calculation of c6 (6) using (6.2). In Table 6.1
the individual terms are listed; note the order of the largest scaled sum-
mands. Yet the sum of the last column is only 1.711 × 10−10 , which
.
when multiplied by 2e6 /π, yields c6 (6) = 0.00000004396.

6.2 A Horner Type Method

This method is based on the simple observation that the series Sr (z)
of (1.2) terminates if z is a non-negative integer. Thus the an (r) can
be recursively determined by successively setting z = 0, 1, 2 . . ., etc.
Specifically, defining
Fr (z) = Γ(z + 1)(z + r + 1/2)−(z+1/2) ez+r+1/2 (2π)−1/2 , (6.3)

74
Chapter 6. Computing the Coefficients

j C2j,12 F6 (j)/e6 C2j,12 F6 (j)/e6

0 1 .810495300735 .810495300735
1 −72 .136735023817 −9.844921714804
2 840 .054363835840 45.665622105199
3 −3584 .029448235581 −105.542476320579
4 6912 .018797607523 129.929063197172
5 −6144 .013277713616 −81.578272459508
6 2048 .010039301705 20.560489891957

Table 6.1: Scaled summands of c6 (6)

then (1.1) may be stated1

1 z z(z − 1)
Fr (z) = a0 (r) + a1 (r) + a2 (r) +···
2 z+1 (z + 1)(z + 2)

from which the coefficients can be efficiently found using a type of

Horner procedure:

a0 = 2Fr (0)

a0 2
a1 = Fr (1) −
2 1

a0 3 4
a2 = Fr (2) − − a1 ,
2 2 1

and in general, for n ≥ 3,

a0 n + 1 n+2 n+3 2n
an (r) = Fr (n) − − a1 − a2 − · · · − an−1 .
2 n n−1 n−2 1
(6.4)
Although not very efficient for computing a single an (r) value, it
proves quite practical and easy to program when a list of coefficients
(for a given value of r) is required, as is most often the case.
1 F (z)is clearly identical to Sr (z). The use of Fr (z) is meant to indicate that
r
this quantity should be computed using the closed form (6.3).

75
Chapter 6. Computing the Coefficients

6.3 Another Recursion Formula

Recall that the Lanczos formula (1.1) is a convergent formula for the
gamma function. As such, it inherits the standard properties of this
function, in particular the fundamental recursion Γ(z + 1) = zΓ(z).
Recall the notation for the rational functions appearing in the infi-
nite series of (1.1): H0 (z) = 1, and for k ≥ 1 an integer,

Γ(z + 1)Γ(z + 1)
Hk (z) =
Γ(z − k + 1)Γ(z + k + 1)
z · · · (z − k + 1)
= . (6.5)
(z + 1) · · · (z + k)

It is a simple exercise to deduce the following relationships:

(N − k)(N + k)
Hk (N − 1) = Hk (N), 0≤k ≤N −1 (6.6)
N2
and
(2N)!
Hk (N) = HN (N), 0≤k≤N . (6.7)
(N − k)!(N + k)!
Now
∞
z+1/2 −(z+r+1/2)
√ X ′
Γ(z + 1) = (z + r + 1/2) e 2π ak (r)Hk (z) , (6.8)
k=0

and
∞
z−1/2 −(z+r−1/2)
√ X ′
zΓ(z) = z(z + r − 1/2) e 2π ak (r)Hk (z − 1) .
k=0

By the fundamental recursion, the right hand sides of these two equa-
tions are equal, from which it follows that
X′
ez(z + r − 1/2)z−1/2 ak (r)Hk (z − 1)
1= X′ . (6.9)
(z + r + 1/2)z+1/2 ak (r)Hk (z)

p
For a0 (r), set z = 0 in equation (6.8) to get a0 (r) = 2e/[π(r + 1/2)] er .
For N ≥ 1, set z = N in (6.9). Both series then terminate, one at the

76
Chapter 6. Computing the Coefficients

N − 1 term while the other at the Nth term. Isolating aN (r) and
simplifying using (6.6) and (6.7) then yields
N −1
" N −1/2 #
X N + r − 1/2 (N − k)(N + k) ak (r)
aN (r) = (2N)! e −1 .
k=0
N + r + 1/2 N(N + r + 1/2) (N − k)!(N + k)!

As with the Horner type of method, this method is not very efficient
for computing a single an (r) value, but is more practical due to recycling
of previous computations when a list of coefficients is required.

6.4 Discrete Fourier Transform

The final method discussed is that of the discrete Fourier transform,
but in the Chebyshev setting. This method takes advantage of the
closed form formula for fE,r (x) determined in Section 4.1.2. The gen-
eral properties of discrete (or finite) Fourier transforms are explained,
followed by the equivalent treatment in the Chebyshev setting.

6.4.1 Finite Fourier Series

The method is based on the principle that the family of exponential
functions {exp (2πikx)}∞ k=−∞ is orthogonal not only with respect to in-
tegration over [0, 1], but also with respect to summation over equispaced
data on this interval. Specifically, in the case of integrals
Z 1
2πinx −2πimx 1 if n = m,
e e dx = (6.10)
0 0 else .

For the discrete analogue, fix an integer N ≥ 2 and consider the values
xj = j/N, j = 1, . . . , N in [0, 1]. Then for 1 ≤ n, m ≤ N,
N
1 X 2πinxj −2πimxj 1 if n = m,
e e = (6.11)
N j=1 0 else .

More generally, for n, m ∈ Z,

N
1 X 2πinxj −2πimxj 1 if n ≡ m (mod N),
e e = (6.12)
N j=1 0 else .

77
Chapter 6. Computing the Coefficients

NowP assume that f (x) has a uniformly convergent Fourier expansion

f (x) = ∞ j=−∞ bj e
2πijx
on [0, 1], and that the values of f are known at
the points xj = j/N, j = 1, . . . , N. Compute approximations to the
Fourier coefficients as
N
X
b̄j = f (xk )e−2πijxk
k=1

and form the discrete Fourier series

N
f¯(x) =
X
b̄j e2πijx .
j=1

Then f¯ = f on the equispaced data points xj , j = 1, . . . , N, and f¯

is a good approximation to f for the points in between. The larger
the number N of sampled data points becomes, the better the approx-
imation. This is clear from equation (6.11) since this expression is an
N-subinterval Riemann sum approximating the integral in (6.10). In
fact, it is possible to derive an explicit relationship between the coeffi-
cients b̄j of the approximation, and the bj of the Fourier series. Since
f¯(xn ) = f (xn ), n = 1, . . . , N, it follows that
N
X ∞
X
2πikxn
b̄k e = bk e2πikxn ,
k=1 k=−∞

again for n = 1, . . . , N. Now multiply both sides by exp (−2πijxn ),

sum over n = 1, . . . , N, and use (6.12) to get
X
b̄j = bk .
k≡j
(mod N)

For this reason, the smoother the function f being approximated, the
more rapidly its Fourier coefficients bj decrease to zero, and hence the
better the one term approximation b̄j ≈ bj .

6.4.2 Finite Chebyshev Series

In the Chebyshev setting, the theory of the previous section translates
as follows: for a given function f (x) on [−1, 1], fix an integer N and

78
Chapter 6. Computing the Coefficients

consider the Chebyshev polynomial TN (x) which has N zeros x1 , . . . , xN

at
π(k − 1/2)
xk = cos .
N
Define the coefficients c̄n by
N
2 X
c̄n = f (xk )Tn (xk )
N k=1

N
2 X π(k − 1/2) nπ(k − 1/2)
= f cos cos .
N k=1 N N

Then the approximation

N
′c̄ T (x)
X
f (x) ≈ n n
n=0

is exact on the N zeros of TN (x). Larger N values produce more accu-

rate estimates of c̄n and hence f (x).
This method can be applied to estimate the coefficients in (1.1)
using the explicit form of fE,r (x) from equation (4.4):
N
2 X
ān = (−1)n c̄n = (−1)n fE,r (xk )T2n (xk ) .
N k=1

This calculation can be economized since fE,r (xk ) and T2n (xk ) are even
with respect to xk , hence, taking N even,
N/2
4 X
n π(k − 1/2) 2nπ(k − 1/2)
ān = (−1) fE,r cos cos .
N k=1 N N

Furthermore, the N values fE,r (x1 ), . . . , fE,r (xN ) need only be com-
puted once since they do not depend on n, the index of the coefficient
being computed.
Selecting N = 20 was sufficient to reproduce a0 , . . . , a5 to six dec-
imals as given in Lanczos’ original paper for r = 1, 1.5, 2, 3. In fact,
N = 20 is really only needed for r = 1, as the error |cn (r) − c̄n | appears

79
Chapter 6. Computing the Coefficients

to drop drastically as r increases. For example, with r = 3 and N = 20,

10 decimal precision is achieved for the first 6 coefficients.
What is not so straightforward, however, is the evaluation of fE,r (x)
to produce the sample function values. The closed form of fE,r (x)
requires evaluation of the Lambert W functions W−1 (y) and W0 (y)
near their branch points y = 0 and y = −1/e. In [5] the authors
discuss the subtleties involved with evaluation near these points. The
numerical experiments noted here were carried out using Maple 8 which
has built-in evaluation routines based on [5].

6.5 Some Sample Coefficients

To sum up this discussion, noted here for the record are coefficients
corresponding to various values of the parameter r. Reproduced in
Table 6.2 are the coefficients given in Lanczos’ original paper. These
values were computed using the Lanczos method.

r=1 r = 1.5 r=2 r=3

a0 /2 +1.4598430249 +2.0844142416 +3.0738046712 +7.0616588080
a1 -0.4606423129 -1.0846349295 -2.1123757377 -6.5993579389
a2 +0.0010544242 +0.0001206982 +0.0386211602 +0.5396522297
a3 -0.0003384921 +0.0001145664 -0.0000510050 -0.0019519669
a4 +0.0001175425 -0.0000176145 +0.0000004776 -0.0000013258
a5 -0.0000506634 +0.0000038119 +0.0000006715 +0.0000002201

Table 6.2: Coefficients as functions of r

6.6 Partial Fraction Decomposition

For the final summary of his results, Lanczos truncates the main for-
mula (1.1) after a finite number of terms and resolves the resulting
rational functions into their constituent partial fractions, resulting in
the formula
" N
#
√ X b k (r)
Γ(z + 1) ≈ (z + r + 1/2)z+1/2 e−(z+r+1/2) 2π b0 (r) + .
k=1
z + k
(6.13)

80
Chapter 6. Computing the Coefficients

This makes for increased efficiency of the formula when multiple eval-
uations are required once a final set of ak (r) coefficients has been com-
puted. For this purpose, it is worthwhile deriving the general formula
for the decomposition.
Recall the infinite series portion of the main formula is

1 z z(z − 1)
Sr (z) = a0 (r) + a1 (r) + a2 (r) +···
2 z+1 (z + 1)(z + 2)
∞
X ′
= ak (r)Hk (z) .
k=0

If this is terminated after the N’thPterm, the goal is to express the

resulting series Sr,N (z) = a0 (r)/2 + N
k=1 ak (r)Hk (z) in the form

N
X bk (r)
Sr,N (z) = b0 (r) + .
k=1
z+k

Observe that for k ≥ 1, bk (r) is the residue of Sr,N (z) at z = −k, while,
upon taking z → ∞,
N
X
b0 (r) = a0 (r)/2 + ak (r) .
k=1

Now
(
(k+j−1)!
(−1)k−j+1 (k−j)![(j−1)!] 2 if 1 ≤ j ≤ k,
Res [Hk (z)]z=−j =
0 else ,

so that
k
X (k + j − 1)! 1
Hk (z) = 1 + (−1)k−j+1 2
.
j=1
(k − j)![(j − 1)!] z + k

To determine bj (r) of equation (6.13), simply sum up the coefficients

of 1/(z + j) in Sr,N (z):
N
X (k + j − 1)!
bj (r) = ak (r)(−1)k−j+1 .
k=j
(k − j)![(j − 1)!]2

81
Chapter 6. Computing the Coefficients

6.7 Matrix Representation

As remarked in Section 6.1, the calculation of the coefficients can be
concisely expressed in matrix form, which, along with the partial frac-
tion decomposition, reduce many of the intermediate calculations to
integer arithmetic thus avoiding some round-off error. Godfrey makes
this observation in [9] and carries through the details. The matrices
required for this purpose are constructed here.
Begin by introducing some notation. Let a(r) be the column vector
of coefficients a0 (r), . . . , aN (r), and let b(r) be the column vector of
partial fraction coefficients b0 (r), . . . , bN (r).
Writing C as the matrix of Chebyshev coefficients C2i,2j , 0 ≤ i, j ≤
N, C has entries

0
 if i < j,
Cij = 1 if i = j = 0,
 (−1)i−j i(i+j)! 4j

(i+j)(2 j)! (i−j)!
else.

Define the partial fraction decomposition matrix to be B, which has

entries 

 0 if i > j,

1/2

if i = j = 0,
Bij =

 1 if i = 0, j > 0,
j−i+1
(−1) (j+i−1)!

else.

 2
(j−i)! ((i−1)!)
as given in Section 6.6.
√
Denote by Fr the column vector of scaled values 2Fr (j)e−r−1/2 ,
0 ≤ j ≤ N, where as in Section 6.1,
Fr (z) = 2−1/2 Γ(z + 1/2)(z + r + 1/2)−z−1/2 exp (z + r + 1/2) .

With this notation, the Fourier coefficients appearing in (1.1) can

be concisely expressed
√ r+1/2
2e
a(r) = CFr
π
while the coefficients for the partial fraction decomposition are simply
√ r+1/2
2e
b(r) = BCFr .
π

82
Chapter 6. Computing the Coefficients

√
The only floating point calculations are in the term ( 2er+1/2 /π)Fr , so
that if multiple sets of coefficients are required for different values of r,
only this term need be recomputed.
Still further efficiency can be realized by observing that the leading
r+1/2
e term of the b(r) coefficients cancels once inserted into (6.13).
Letting u(z) denote the vector of rational functions 1, (z+1)−1 , . . . , (z+
1)−k , equation (6.13) becomes simply
r z+1/2
e z + r + 1/2
Γ(z + 1) ≈ 2 uT (z)BCFr . (6.14)
π e

83
Chapter 7
Lanczos’ Limit Formula

The final result of Lanczos’ paper [14] is the beautiful limit formula,
stated here as

Theorem 7.1. For Re(z) ≥ 0,

∞
X ′ 2 /r
z
Γ(z + 1) = 2 lim r (−1)k e−k Hk (z) . (7.1)
r→∞
k=0

This result1 seems to be a departure from the practical nature of the

rest of the paper, and the connection to the previous formula (1.1) is
not an obvious one. Indeed, Lanczos’ reviewer [29] notes “Exceedingly
curious is the additional remark, made by the author without any proof
. . . ”. Lanczos offers little insight into how this result follows from
his main derivation, other than the two sentences which precede the
formula:

If r grows to infinity, we obtain a representation of the facto-

rial function which holds everywhere in the complex plane.
In this case we are able to give the coefficients of the se-
ries . . . in explicit form, due to the extreme nature of the
function v r .

In this chapter the validity of (7.1) is proved in detail for Re(z) ≥ 0.

1 Note that the formula is misstated in the review [29]. There the term r z is
incorrectly stated as r 2 .

84
Chapter 7. Lanczos’ Limit Formula

7.1 Some Motivation

As noted, it is difficult to say what lead Lanczos to state in a matter
of fact fashion that (7.1) should be true. Some insight can be gained,
however, by examining the behaviour of the integral in equation (3.6)
for large positive r. To this end, write (3.6) as
Z π/2
z+1/2 −(z+r+1/2) 2v r sin θ
Γ(z + 1/2) = (z + r + 1/2) e cos2z θ dθ
−π/2 log v

z+1/2 Z π/2
2 (v/e)r sin θ

z + r + 1/2
= cos2z θ dθ .
e −π/2 log v
A rescaling of r to er produces the r z term of equation (7.1) and elim-
inates the exponential term in the limit, thus
z+1/2 Z π/2
r 1/2 (v/e)er sin θ

z z + er + 1/2
Γ(z + 1/2) = 2r cos2z θ dθ
er −π/2 log v

z+1/2 Z π/2
z z + er + 1/2
= 2r cos2z θ r 1/2 gr (θ) dθ say,
er −π/2
(7.2)

where
sin θ
gr (θ) = (v/e)er.
log v
Now gr (θ) converges point-wise to zero rapidly on [−π/2, π/2) with
increasing r, but equals one at θ = π/2. As a result, as r increases,
(r/π)1/2 gr (θ) takes on the appearance of a the Gaussian distribution
(r/π)1/2 exp [−r(θ − π/2)2]. Plots of exp [−r(θ − π/2)2 ] − gr (θ) for in-
creasing values of r show a convincing fit; see Figure 7.1.
Furthermore,
Z ∞
1 2 2
√ cos (2kθ)r 1/2 e−r(θ−π/2) dt = (−1)k e−k /r .
π −∞
so that if gr (θ) is replaced by exp [−r(θ − π/2)2 ] and the bounds of
integration [−π/2, π/2] are extended to (−∞, ∞), the integral gives
2
rise to the terms e−k /r upon writing cos2z θ as its Fourier series.

85
Chapter 7. Lanczos’ Limit Formula

0.08

r=1

r=6

r=10

− π2 0 π
2

Figure 7.1: exp [−r(θ − π/2)2 ] − gr (θ) for r = 1, 6, 10

These ideas lead to the following outline of the proof of equa-

tion (7.1). For large r and fixed z,
z+1/2 Z π/2
z z + er + 1/2
Γ(z + 1/2) = 2r cos2z θ r 1/2 gr (θ) dθ
er −π/2

Z π/2
z 2
≈ 2r cos2z θ r 1/2 e−r(θ−π/2) dθ
−π/2

Z ∞
z 2
≈r cos2z θ r 1/2 e−r(θ−π/2) dθ
−∞

86
Chapter 7. Lanczos’ Limit Formula

∞
!
Z ∞ ′ 2 Γ(z + 1/2) 2
X
= rz √ Hk (z) cos (2kθ) r 1/2 e−r(θ−π/2) dθ
−∞ 0
π Γ(z + 1)

∞
z Γ(z + 1/2) X′ 2
= 2r (−1)k e−k /r Hk (z) ,
Γ(z + 1) 0

where the ≈ symbol becomes equality as r → ∞. The remainder of the

chapter is dedicated to making precise the details of this argument.

7.2 Proof of the Limit Formula

The steps of the proof will be justified by a series of smaller lemmas.
Begin with a definition:

Definition 7.1. For r > 0, define

1/4
1
δr = . (7.3)
r log r

Clearly δr → 0 as r → ∞.

It will also be convenient to introduce some notation for the various

intervals of integration:

Definition 7.2. For δr < π, define the intervals

hπ πi
Ir = − δr , ,
2 2
h π π i
Jr = − , − δr , and
2 2
π i
Kr = −∞, − δr .
2

The first lemma is straightforward:

87
Chapter 7. Lanczos’ Limit Formula

Lemma 7.1. The function

(
0 if θ = 0,
h(θ) = sin θ
log v
else

is continuous on [−π/2, π/2], and is thus bounded.

Proof of Lemma 7.1: On (−π/2, π/2], h(θ) = fer (θ)v −er 2−1/2 , which
is continuous by Theorem 4.3. At θ = −π/2,

lim h(θ) = 0 = h(−π/2) .

θ→−π/2

The next lemma shows that the integral in equation (7.2) is con-
centrated on the small interval Ir containing π/2:

Lemma 7.2. Suppose Re(z) ≥ 0, and let

Z v er sin θ
I1 (r, z) = cos2z θ r 1/2 dθ .
Jr e log v

Then
lim r z I1 (r, z) = 0 .
r→∞

Proof of Lemma 7.2: First, recall the calculus result

sin2 x
lim =1,
x→0 x2

so that for |x| sufficiently small,

1
sin2 x ≥ x2 . (7.4)
2

88
Chapter 7. Lanczos’ Limit Formula

Now bound I1 (r, z) as follows:

Z v er sin θ
z 2z 1/2
r cos θ r dθ

Jr e log v

Z v er
z 2z 1/2
≤ C1 r
cos θ r dθ for some C1 > 0 by Lemma 7.1
Jr e

er
σ+1/2 v(π/2 − δr )
≤ C1 r |Jr | since v is increasing on Jr
e

≤ C2 r σ+1/2 eer[log v(π/2−δr )−1]

2
= C2 r σ+1/2 e−er cos (π/2−δr )/v
by the implicit definition of v(θ)

2
≤ C2 r σ+1/2 e−er cos (π/2−δr )/e

2
= C2 r σ+1/2 e−r sin δr

2
≤ C2 r σ+1/2 e−rδr /2 for δr sufficiently small, by (7.4)

√
= C2 r σ+1/2 e−(1/2) r/ log r

−→ 0 as r → ∞

The next lemma justifies the replacement of gr (θ) by exp [−r(θ − π/2)2 ]
on Ir :

Lemma 7.3. Suppose Re(z) ≥ 0, and let

v er sin θ
Z
2z 1/2 −r(θ−π/2)2
I2 (r, z) = cos θ r −e dθ .
Ir e log v

89
Chapter 7. Lanczos’ Limit Formula

Then
lim r z I2 (r, z) = 0 .
r→∞

Proof of Lemma 7.3: Begin with

Z
z
r 2z 1/2 v er sin θ −r(θ−π/2)2

cos θ r −e dθ

Ir e log v

v er sin θ r(θ−π/2)2
Z
σ 2σ 1/2 −r(θ−π/2)2

≤r cos θ r e
e e − 1 dθ .
Ir log v

Letting χr denote the indicator function on Ir , it is enough to show

that as r → ∞

(i)
v er sin θ r(θ−π/2)2

e e − 1 χr −→ 0
log v
uniformly, and

(ii) that Z
σ 2
r cos2σ θ r 1/2 e−r(θ−π/2) dθ
Ir

remains bounded.

For (i), since limθ→π/2− sin θ/ log v = 1, it is enough to show that

h v er i h i
r(θ−π/2)2 er(log v−1)+r(θ−π/2)2
e − 1 χr = e − 1 χr −→ 0
e
uniformly, which boils down to showing that

er(log v − 1) + r(θ − π/2)2 χr −→ 0

uniformly. Expanding e(log v − 1) + (θ − π/2)2 about θ = π/2 gives

e(log v − 1) + (θ − π/2)2

1 1 4 1 135 60
(θ − π/2)6 + O (θ − π/2)6 .

= − (θ − π/2) + −4 − 2 +
3 e 90 e e

90
Chapter 7. Lanczos’ Limit Formula

Now the coefficient (1/3 − 1/e) < 0, so that for |θ − π/2| sufficiently
small,
−(θ − π/2)4 ≤ e(log v − 1) + (θ − π/2)2 ≤ 0 ,
so that, for r sufficiently large,

−(θ − π/2)4 ≤ e(log v − 1) + (θ − π/2)2 χr ≤ 0 ,

whence

−r(θ − π/2)4 ≤ er(log v − 1) + r(θ − π/2)2 χr ≤ 0 .

The term (θ − π/2) is at most δr in magnitude on Ir , so

−rδr4 ≤ er(log v − 1) + r(θ − π/2)2 χr ≤ 0 ,

and hence

1
≤ er(log v − 1) + r(θ − π/2)2 χr ≤ 0 .

−r
r log r

Now let r → ∞ to get

er(log v − 1) + r(θ − π/2)2 χr −→ 0 uniforly.

This proves (i).

For (ii), first observe that about θ = π/2,

cos2σ θ ≤ (θ − π/2)2σ .

Thus
Z
σ 2
r cos2σ θ r 1/2 e−r(θ−π/2) dθ
Ir

Z
σ 2
≤r (θ − π/2)2σ r 1/2 e−r(θ−π/2) dθ
Ir

Z ∞
σ 2
≤r (θ − π/2)2σ r 1/2 e−r(θ−π/2) dθ .
π/2

91
Chapter 7. Lanczos’ Limit Formula

Now make the change of variables u = r(θ − π/2)2 , du = 2r(θ − π/2)dθ.

This gives
Z ∞
σ 2
r (θ − π/2)2σ r 1/2 e−r(θ−π/2) dθ .
π/2

∞
1 u σ−1/2
Z
= r σ−1 r 1/2 e−u du
2 0 r

∞
1
Z
= uσ−1/2 e−u du
2 0

1
= Γ(σ + 1/2) .
2

That is,
1
Z
σ 2
r cos2σ θ r 1/2 e−r(θ−π/2) dθ ≤ Γ(σ + 1/2) .
Ir 2

This proves (ii) and completes the proof.

Lemma 7.4. Suppose Re(z) ≥ 0, and let

Z
2
I3 (r, z) = cos2z θ r 1/2 e−r(θ−π/2) dθ .
Kr

Then
lim r z I3 (r, z) = 0 .
r→∞

92
Chapter 7. Lanczos’ Limit Formula

Proof of Lemma 7.4:

Z
1/2 −r(θ−π/2)2
z 2z

r cos θ r e dθ

Kr

Z
σ 2
≤r r 1/2 e−r(θ−π/2) dθ
Kr

Z ∞
σ 2
=r √
e−u du ,
rδr

√
this last line a result of the change of variables u = r(θ − π/2). This
can be expressed in terms of the error function
Z x
2 2
erf(x) = √ e−t dt
π 0
2
e−x
≈ 1− √ as x → ∞ .
πx

Picking up where we left off then we have

Z ∞ √
σ −u2 σ π
√
r √ e du = r 1 − erf( rδr )
rδr 2

√
−( rδr )2
σ1e
≈r √
2 rδr

1 1/2
= r σ−1/4 (log r)1/4 e−(r/ log r)
2

−→ 0 as r → ∞ ,

completing the proof.

93
Chapter 7. Lanczos’ Limit Formula

Lemma 7.5. For integer k,

Z ∞
1 2 2
√ cos (2kθ)r 1/2 e−r(θ−π/2) dθ = (−1)k e−k /r . (7.5)
π −∞

Proof of Lemma 7.5: The integral in (7.5) is the real part of

Z ∞
1 2
√ r 1/2 ei2kθ−r(θ−π/2) dθ
π −∞

∞
1
Z
i2k+πr 2
k −k 2 /r
= (−1) e √ r 1/2 e−r(θ− 2r
)
dθ upon completing the square
π −∞

∞
1
Z
k −k 2 /r 2
= (−1) e √ r 1/2 e−rθ dθ with a simple contour shift
π −∞

2 /r
= (−1)k e−k .

We are now in a position to prove (7.1).

Proof of Theorem 7.1: Recall equation (7.2) which is true for all
r ≥ 0:
z+1/2 Z π/2
r 1/2 (v/e)er sin θ

z z + er + 1/2
Γ(z+1/2) = 2r cos2z θ dθ .
er −π/2 log v
From the definitions of I1 (r, z), I2 (r, z) and I3 (r, z) in Lemmas 7.2, 7.3
and 7.4, respectively,
z+1/2 Z π/2
r 1/2 (v/e)er sin θ

z z + er + 1/2
2r cos2z θ dθ
er −π/2 log v

z+1/2
z z + er + 1/2
= 2r
er
1 ∞
Z
2z 1/2 −r(θ−π/2)2
× I1 (r, z) + I2 (r, z) − 2I3 (r, z) + cos θ r e dθ .
2 −∞

94
Chapter 7. Lanczos’ Limit Formula

Now let r → ∞ and apply Lemmas 7.2, 7.3 and 7.4 to get
Z ∞
z1 2
Γ(z + 1/2) = 2 lim r cos2z θ r 1/2 e−r(θ−π/2) dθ .
r→∞ 2 −∞
Finally, replace cos2z θ with its Fourier series (see Section 3.3)
∞
X ′ 2 Γ(z + 1/2)
2z
cos θ = √ Hk (z) cos (2kθ)
0
π Γ(z + 1)
and integrate term by term using Lemma 7.5 to get
∞
z Γ(z+ 1/2) X′ 2
Γ(z + 1/2) = 2 lim r (−1)k e−k /r Hk (z) .
r→∞ Γ(z + 1) 0

Canceling the Γ(z + 1/2) terms and moving Γ(z + 1) to the left hand
side completes the proof.

7.3 Additional Remarks

In his original paper [14], Lanczos claims that the limit in Theorem 7.1
converges for all z away from the negative integers and defines the
gamma function in this region. The proof given here, however, does
not extend to the left hand plane.
From his comments following the statement of the theorem, he
seems to base the limit formula on the limiting behaviour of the individ-
ual coefficients ak (r), not on the limiting behaviour of the integral (7.2).
Specifically, the comparison appears to be
∞
√ z+1/2 −(z+r+1/2)
X ′
Γ(z + 1) = 2π (z + r + 1/2) e ak (r)Hk (z)
k=0

∞ r
z+1/2 X
z z + er + 1/2 ′ πr ak (er)
= 2r Hk (z)
er 2 eer
k=0

∞
X ′ 2 /r
z
≈ 2r (−1)k e−k Hk (z) .
0

95
Chapter 7. Lanczos’ Limit Formula

This suggests that for large r,

r
πr ak (er) 2

er
∼ (−1)k e−k /r ,
2 e
or, by rescaling r,
r
r−ek 2 /r 2e
Q(k, r) = (−1)k e (ak (r))−1 ∼ 1 .
πr
Plots of log Q(k, r) for r = 10, 20 and 50 support this conjecture; refer
to Figure 7.2, 7.3 and 7.4, taking note of the different vertical scales.
90

−90
0 k 24

Figure 7.2: log Q(k, 10), k = 0, . . . , 24

96
Chapter 7. Lanczos’ Limit Formula

35

−35
0 k 24

Figure 7.3: log Q(k, 20), k = 0, . . . , 24

−2
0 k 24

Figure 7.4: log Q(k, 50), k = 0, . . . , 24

97
Chapter 8
Implementation & Error
Discussion

In order to use the main formula (1.1) in practice to compute values

of Γ(z + 1), the series P (1.2) is truncated and an estimate of the tail
∞
of the series ǫr,n (z) = k=n+1 ak (r)Hk (z) is therefore required. The
error function ǫr,n (z) is essentially the relative error, and in this chapter
uniform bounds on this function are established for z in the right-half
complex plane.
The chapter begins with a brief discussion of the various notions
of error which arise in approximation theory. Following the introduc-
tion of some new notation and an important theorem, Lanczos’ own
error estimates and those of Luke [17] are examined. The final sec-
tion deals with numerical investigations motivated by the observation
that the error function ǫr,n (z) is the tail of an absolutely and uniformly
convergent series, and as such, possesses the desirable properties of
an analytic function. In particular, the maximum modulus of the error
function in the right half plane will be shown to occur on the imaginary
axis. Furthermore, the value of this maximum can be easily estimated
empirically with the mapping t → it/(1 − t) of the interval [0, 1) onto
the positive imaginary axis. Finally, the limiting behaviour of ǫr,n (z)
as z → ∞ is examined which leads to error bounds much improved on
those given by Lanczos in [14].
The principal results of this chapter are the proof that the maxi-
mum of |ǫr,n (z)| in the right half plane occurs on the imaginary axis
(Theorem 8.1), and the optimal formulas resulting from the method of
empirically bounding ǫr,n (z) given in Section 8.5.

98
Chapter 8. Implementation & Error Discussion

8.1 A Word or Two about Errors

Before we begin, we should clarify the various notions of “error” present
in the literature. The definition and use of absolute error is a universal
standard, and the distinction between absolute and relative error is
fairly clear. The term “relative error” itself, however, is given slightly
different treatment depending on the author. At the end of the day,
these different treatments amount to essentially the same thing, but the
subtle differences between them are worth noting and so are pointed
out here.

8.1.1 Absolute and Relative Error

This discussion will be in the context of the gamma function. As such,
let G(1 + z) denote an approximation to Γ(1 + z), and let ǫα denote
the difference between the approximation and the actual value of the
function:
ǫα (z) = Γ(1 + z) − G(1 + z) .
This is often termed the “absolute error”, even without taking the
modulus of ǫα . The relative error is here denoted ǫρ (z), and it is defined
according to the commonly accepted notion of (actual-estimate)/actual:

Γ(1 + z) − G(1 + z)
ǫρ (z) =
Γ(1 + z)

ǫα (z)
= . (8.1)
Γ(1 + z)

From the definition of ǫρ (z), and assuming |ǫρ (z)| < 1, Γ(1 + z) may
be expressed

G(1 + z)
Γ(1 + z) = (8.2)
1 − ǫρ (z)

= G(1 + z)[1 + ǫρ (z) + ǫρ (z)2 + · · · ] (8.3)

≈ G(1 + z)(1 + ǫρ (z)) . (8.4)

99
Chapter 8. Implementation & Error Discussion

It is common in the literature to encounter expressions of the type

Γ(1 + z) = G(1 + z)(1 + ǫ(z)), where ǫ(z) is termed the relative error.
This is not, strictly speaking, relative error in the sense of (8.1), but
comparison with equations (8.3) and (8.4) shows that it is essentially
the same (up to first order).
One may ask why relative error is of concern in the first place, when
it seems more natural to ask only that estimated quantities be within
a prescribed accuracy of their true values. The answer is that relative
error is a measure of error as a proportion of the true answer, so that nu-
merically, it reports the number of reliable digits in our approximation.
This is an especially important consideration when performing floating
point calculations in which the maximum number of digits is fixed. To
appreciate the distinction, consider an example: suppose two different
approximations x1 = 123 and x2 = 0.123 are computed, both with a
relative error of ǫρ = 0.1. Then the actual errors are xk ǫρ /(1 − ǫρ ),
k = 1, 2, which are approximately 14 and 0.014, respectively. This
means that for both x1 and x2 , the first non-zero digit is good while
the second may be in question. By contrast, if x1 and x2 both have
absolute error 0.1, then the estimate for x2 = 123 appears acceptable,
while that for x2 = 0.123 is clearly a poor one.

8.1.2 Lanczos’ Relative Error

The Lanczos notion of relative error is quite different from both (8.1)
and (8.4). In his work, he writes
√
Γ(z + 1) = (z + r + 1/2)z+1/2 e−(z+r+1/2) 2π [Sr,n (z) + ǫr,n (z)] (8.5)

where Sr,n (z) is the series (1.2) truncated after a finite number of terms,
and he calls ǫr,n (z) the relative error. It turns out that Sr,n (z) ≈ 1 in
the right half plane, so that Lanczos’ notion is very close to (8.4). The
relationship between ǫr,n (z) and (8.1) is
ǫr,n (z)
ǫρ = .
Γ(z + 1)(z + r + 1/2)−(z+1/2) ez+r+1/2 (2π)−1/2
Spouge makes this same observation in [27] and shows that for Re(z) ≥
0, the denominator is bounded below by
e 1/2
|Γ(z + 1)(z + r + 1/2)−(z+1/2) ez+r+1/2 (2π)−1/2 | ≥ ,
π

100
Chapter 8. Implementation & Error Discussion

where upon selecting r = z = 0 the bound in this inequality is sharp.

Thus π 1/2
|ǫρ | ≤ |ǫr,n (z)| , (8.6)
e
p .
and since (π/e) = 1.075, bounding relative error in the Lanczos sense
gives essentially the same bound for ǫρ , the standard notion of relative
error. For the purposes of this study z will be restricted to Re(z) ≥ 0
and as such the Lanczos notion of relative error will be used.
It should be pointed out, however, that with ǫr,n (z) something is
lost if z is moved to the left of the imaginary axis. As we saw in
Chapter 5, the formula (1.1) extends to the region Re(z) > −r not
including the negative integers. Lanczos’ relative error ǫr,n (z) has poles
at the negative integers in this case, and becomes unbounded. ǫρ (z), on
the other hand, is analytic in this region, and it is therefore possible to
give meaning to the relative error in the approximation near the poles.

8.1.3 Restriction to Re(z) ≥ 0

To conclude the introductory remarks about errors, the following com-
mon practice specific to the gamma function is important. Recall the
reflection formula (2.10) which permits the evaluation of gamma in the
left half plane once its values in Re(z) ≥ 0 are known:
πz
Γ(1 + z)Γ(1 − z) = .
sin πz
Already this says that methods used to estimate Γ(1 + z) need only be
valid on Re(z) ≥ 0.
Now suppose Γ(1 + z) is computed with a relative error of ǫρ , that
is,
G(1 + z)
Γ(1 + z) = .
1 − ǫρ
Then by the reflection formula,
πz 1 − ǫρ
Γ(1 − z) =
sin πz G(1 + z)
πz 1 1
=
sin πz G(1 + z) 1 + ǫρ + ǫ2ρ + · · ·
πz 1 1
≈ .
sin πz G(1 + z) 1 + ǫρ

101
Chapter 8. Implementation & Error Discussion

In other words, to compute Γ(1 − z) to a relative error of −ǫρ , it is

sufficient to compute Γ(1 + z) to a relative error of ǫρ , and vice versa.
Thus when using a particular method to compute Γ(1+z), it is sufficient
to be able to bound relative error for Re(z) ≥ 0.

8.2 Some Notation and a Theorem

This section lists notation which will be used throughout the remainder
of this study, and concludes with the proof of an important theorem.
Some of the functions noted here have been introduced previously, but
are repeated here for easy reference.

1. Recall the definition of Hk (z):

H0 (z) = 1 ,

and for k ≥ 1,
Γ(z + 1)Γ(z + 1)
Hk (z) =
Γ(z − k + 1)Γ(z + k + 1)

z · · · (z − k + 1)
=
(z + 1) · · · (z + k)

Observe that Hk (z) is a meromorphic function with simple poles

at z = −1, . . . , −k and simple zeros at z = 0, . . . , k − 1. Also,
Hk (z) → 1 as |z| → ∞.

2. Write (1.1) as
√
Γ(z + 1) = 2π (z + r + 1/2)z+1/2 e−(z+r+1/2) [Sr,n (z) + ǫr,n (z)]

where n
X ′
Sr,n (z) = ak (r)Hk (z)
k=0

and ∞
X
ǫr,n (z) = ak (r)Hk (z) .
k=n+1

102
Chapter 8. Implementation & Error Discussion

3. Let

Fr (z) = Γ(z + 1)(z + r + 1/2)−(z+1/2) ez+r+1/2 (2π)−1/2 .

By Stirling’s formula, Fr (z) → 1 as |z| → ∞ in any half plane

Re(z) ≥ σ.
Note that
ǫr,n (z) = Fr (z) − Sr,n (z), (8.7)
and that n
X ′
lim ǫr,n (z) = 1 − ak (r) ,
|z|→∞
k=0

again in any half plane Re(z) ≥ σ, in particular for z in the right

half plane Re(z) ≥ 0. In view of this last limit, the error at
infinity is denoted
n
X ′
ǫ∞
r,n =1− ak (r) .
k=0

4. Let n
X ′
ηr,n (θ) = fE,r (θ) − ak (r) cos (2kθ) ,
k=0

the error in the (n + 1)-term Fourier series approximation of

fE,r (θ), where fE,r (θ) has the meaning of (3.7). Then ηr,n (−π/2)
is the Fourier error at θ = −π/2, which is, incidentally, equal to
ǫr,n (−1/2).

From equation (8.7) we state

Lemma 8.1. For r ≥ 0, ǫr,n (z) is analytic as a function of z in Re(z) >

−1/2.

Proof of Lemma 8.1: In this region Re(z + r + 1/2) > 0 so that

Fr (z) is analytic there, as are the functions Hk (z), k ≥ 0, and hence by
equation (8.7), so is ǫr,n (z).

The following lemma is a consequence of the maximum modulus

principle:

103
Chapter 8. Implementation & Error Discussion

Lemma 8.2. Suppose f is analytic in the right half plane Ω = {Re(z) ≥

σ} and has the property

lim |f (z)| = C < ∞ . (8.8)

|z|→∞
z∈Ω

Let B = supz∈Ω |f (z)|. Then B = supRe(z)=σ |f (z)|.

Proof of Lemma 8.2: Since f is analytic on Ω, by (8.8) it is bounded

so B exists, and we have C ≤ B. Now there are two possibilities:

(i) If B = C, then B = limt→∞ |f (σ + it)|; otherwise

(ii) B > C. For this second case, for any compact K ⊂ Ω, denote
by MK the maximum of |f | on K, and consider the sequence of
closed right-half semi-disks {Ωj }∞
j=1 centered at the origin, each
Ωj having radius j and diameter on the line Re(z) = σ. By the
maximum modulus principle,

MΩj = M∂Ωj ,

and M∂Ωj increases monotonically to B. Write the boundary of

Ωj as the union of two closed segments

∂Ωj = Ij ∪ Aj ,

where Ij is the diameter and Aj the semicircular arc of ∂Ωj . Now

lim M∂Ωj = B ,
j→∞

so that
lim max {MIj , MAj } = B ,
j→∞

but
lim MAj = C < B .
j→∞

Thus
lim MIj = B
j→∞

completing the proof.

104
Chapter 8. Implementation & Error Discussion

Theorem 8.1. For r ≥ 0 and Re(z) ≥ 0, the supremum of |ǫr,n (z)| in

this region occurs on the line z = it, possibly at z = i∞.

Proof of Theorem 8.1: Recall that

n
X ′
lim ǫr,n (z) = 1 − ak (r) .
|z|→∞
k=0

ǫr,n (z) is analytic by Lemma 8.1; now apply Lemma 8.2.

Theorem 8.1 is significant since it reduces the problem of bounding

|ǫr,n (z)| in the right half plane to that of bounding it on the imaginary
axis.

8.3 Lanczos’ Error Estimates

The following passage taken from [14] is the only commentary Lanczos
makes concerning errors in the approximating formula:

A good check on the accuracy of the truncated Fourier

series is provided by evaluating the approximate value of
fr (θ)1 at θ = −π/2, that is by forming the alternate sum
1 π er
ρ0 − ρ1 + ρ2 − ρ3 + · · · = fr − = √ .
2 2 2
We know from the theory of the Fourier series that the max-
imum local error (after reaching the asymptotic stage) can
be expected near to the point of singularity. Let this error
be η. Then a simple estimation shows that the influence
of this error on the integral transform (16) (for values of
z which stay within the right complex half plane), cannot
be greater than (π/2)η. Thus we can give a definite error
bound for the approximation obtained.
1 Inhis paper, Lanczos uses γ instead of r to denote the free parameter ap-
pearing in the formula. The variable r is used here to avoid confusion with
Euler’s constant.

105
Chapter 8. Implementation & Error Discussion

The “integral transform (16)” referred to here is the integral in (3.6):

Z +π/2
cos2z (θ)fr (θ) dθ ,
−π/2

and Lanczos’ comments relate the error in truncating the series in (1.1)
after the nth term, ǫr,n (z), to the error ηr,n (θ) in the (n+1)-term Fourier
series approximation of fE,r (θ). Based apparently on this observation,
he gives uniform error bounds for Re(z) ≥ 0 as listed in Table 8.1.

n r |ǫr,n (z)| <

1 1 0.001
1 1.5 0.00024
2 2 5.1 × 10−5
3 2 1.5 × 10−6
3 3 1.4 × 10−6
4 4 5 × 10−8
6 5 2 × 10−10

Table 8.1: Lanczos’ uniform error bounds

The steps connecting Lanczos’ observation about ηr,n (θ) to the uni-
form error bounds in the table elude me, and it is unclear how he makes
the leap from one to the other. One of his error bounds does appear
incorrect, although this may simply be a matter of round-off. For r = 4
and n = 4, he gives |ǫ4,4 (z)| < 5×10−8 . It is a simple matter to compute
directly
4
a0 (4) X
lim |ǫ4,4 (it)| = 1 − − ak (4)
t→∞ 2 k=1
.
= 5.3 × 10−8 .

His assertion about maximum errors in Fourier approximation oc-

curring “near to the point of singularity” is questionable as well, if not
incorrect altogether. It seems to rely on some sort of Gibbs effect prin-
ciple, and would be accurate if fE,r (θ) had a finite jump discontinuity at
θ = −π/2. In this case, however, fE,r (θ) is continuous at the endpoints

106
Chapter 8. Implementation & Error Discussion

n r θr,n
1 1 −π/2
1 1.5 0
2 2 −π/2
3 2 ≈ −0.72
3 3 −π/2
4 4 0
6 5 −π/2

Table 8.2: Location of maximum of ηr,n (θ)

n r |ǫ∞
r,n |
1 1 8.0 × 10−4
1 1.5 2.2 × 10−4
2 2 5.0 × 10−5
3 2 9.1 × 10−7
3 3 1.1 × 10−6
4 4 5.3 × 10−8
6 5 1.9 × 10−10

Table 8.3: Errors at infinity

and has an infinite jump discontinuity in one of its higher derivatives

(its ⌊2r⌋ + 1 derivative to be precise). Numerical checks show that
Lanczos’ assertion only holds in four of the seven cases noted. Specif-
ically, letting θr,n denote the least value of θ in [−π/2, π/2] where the
maximum of ηr,n (θ) occurs, we find (empirically) the values listed in Ta-
ble 8.2. Indeed, it is not difficult to construct examples in which ηr,n (θ)
is zero at θ = −π/2. Take for example n = 1 and r = 1.500773 · · · .
In this case the maximum of |ηr,n (θ)| is approximately 0.00024 and it
occurs at θ = 0.
One may surmise that Lanczos’ bounds are, for each choice of r and
n, simply the error at infinity, ǫ∞
r,n . A little calculation shows this is
close, but not the case in all instances; see Table 8.3. The apparent
coincidence is likely due to the fact that for many choices of r, the
maximum of |ǫr,n (z)| does occur as |z| → ∞. This is not the case in
general though, as will soon be seen.

107
Chapter 8. Implementation & Error Discussion

Despite the lack of explanation, Lanczos’ uniform error bounds are

all correct, except for the slight difference in the r = n = 4 case as
noted. Empirical checks show that his bounds are not tight, though
there is little room for improvement.

8.4 Luke’s Error Estimates

Except for verbatim quotes of Lanczos’ own estimates [7] [21], the only
other work in this direction to be found in the literature is that of Luke
in [17, p.31]. There the author makes observations similar to Lanczos’
concerning the maximum of ηr,n (θ) near singularities of fE,r (θ), and he
uses this notion to bound ǫr,n (z) for σ = Re(z) > 0. Luke’s bound is
reproduced here.
From the derivation in Section 3,
∞
X
ǫr,n (z) = ak (r)Hk (z)
k=n+1
∞ Z π/2
1 Γ(z + 1) X
= √ ak (r) cos2z θ cos 2kθ dθ
π Γ(z + 1/2) k=n+1 −π/2

π/2
1 Γ(z + 1)
Z
= √ cos2z θηr,n (θ) dθ
π Γ(z + 1/2) −π/2

Thus, assuming |ηr,n (θ)| is maximized at θ = −π/2,

1 Γ(z + 1) π/2
Z
|ǫr,n (z)| ≤ √ cos2σ θ|ηr,n (−π/2)| dθ
π Γ(z + 1/2)
−π/2

1 ηr,n (−π/2)Γ(z + 1) √ Γ(σ + 1)Γ(σ + 1/2)

= √ π Γ(σ + 1)Γ(σ + 1)
π Γ(z + 1/2)

ηr,n (−π/2)Γ(z + 1) Γ(σ + 1/2)
= Γ(σ + 1) .

Γ(z + 1/2)

Unfortunately, this argument breaks down when one considers that

ηr,n (−π/2) has zeros as a function of r, as was pointed out in the
discussion of the Lanczos’ error bounds. Otherwise, we would find

108
Chapter 8. Implementation & Error Discussion

ourselves in the most agreeable situation of having a closed form formula

for Γ(z+1) on the right half plane. Furthermore, this error bound is not
uniform since it grows without bound as |z| → ∞ along lines parallel
to the imaginary axis.
This subtle role of r on the error bounds motivates a closer exami-
nation of ǫr,n (z) as a function of r.

8.5 Numerical Investigations & the Zeros

of ǫ∞
r,n

When terminating the series (1.2) at a finite number of terms, other

authors [7] [27] suggest an optimal strategy of choosing the truncation
order n based on the value of the parameter r. The guiding principle
in this approach seems to be that large r values produce small relative
errors. However, this approach does not take full advantage of many of
the nice properties of the error as a complex analytic function. Instead,
approach the problem the other way around: for a given truncation or-
der n, what value of r minimizes the uniform error bound? This section
presents the results of a series of numerical investigations undertaken
to shed some light on the the dependence of the relative error as a
function of r.

8.5.1 The Motivating Idea

The idea to examine ǫr,n (z) as a function of r stems from the following
line of reasoning: by Theorem 8.1, for a fixed r and n, the maximum
MΩ of |ǫr,n (z)| in the right half plane occurs on the imaginary axis,
possibly at z = ±i∞. In fact, by the Schwartz reflection principle, one
need only consider the positive imaginary axis. Using Stirling’s series
as a guide, if possible, select r = r∗ so that the corresponding error at
infinity, ǫ∞
r∗ ,n , is zero. This forces the location of the maximum onto a
bounded segment of the imaginary axis, hopefully close to the origin.
This gives two possibilities (for this r = r∗ ):

109
Chapter 8. Implementation & Error Discussion

case (i): MΩ = 0

This case is impossible, for otherwise the result would be a closed form
expression
√
Γ(z + 1) = 2π (z + r∗ + 1/2)z+1/2 e−(z+r∗ +1/2) Sr∗ ,n (z)
which could be continued analytically to the slit plane C \ (−∞, −r∗ −
1/2), but which would remain bounded near −(n + 1), −(n + 2), . . .,
contrary to the unboundedness of Γ at the negative integers.

case (ii): MΩ > 0

In this case the maximum of |ǫr∗ ,n (z)| occurs on a bounded segment of

the imaginary axis; the question remains: how to locate it? Experi-
mentally, for each n = 0, . . . , 10, zeros r∗ were found and |ǫr∗ ,n (z)| on
the segment [0, 2in] was examined. In each case the error function was
found to have a single local maximum on this set. Could there be other
local maxima further up the imaginary axis? The difficulty is that
|ǫr∗ ,n (z)| in the form (8.7) is a complicated function, and maximizing it
on the positive imaginary axis using analytical methods is not easy.
The problem is partially overcome with the following observation:
under the transformation z = it/(1 − t), 0 ≤ t < 1, which maps the
unit interval onto the positive imaginary axis, the rational functions
Hk (z) take the form
k−1
Y z(t) − j
Hk (z(t)) =
j=0
z(t) + j + 1

k−1
Y it/(1 − t) − j
=
j=0
it/(1 − t) + j + 1

k−1
Y t(i + j) − j
= .
j=0
t(i − j − 1) + j + 1

The Hk (z(t)) expressions in this last form are numerically stable and
readily computable for 0 ≤ t < 1. The upshot of this observation is that

110
Chapter 8. Implementation & Error Discussion

under the assumption that the first few terms of series ǫr,n (z) contain
the bulk of the error, which turned out to be the case in practice, the
error can be easily estimated as
M
X
ǫr∗ ,n (z(t)) ≈ ak (r∗ )Hk (z(t))
k=n+1

for M not too large. The problem of bounding |ǫr∗ ,n (z)| in the right half
plane reduces to that of estimating the maximum of M
P
k=n+1 ak (r∗ )Hk (z(t))
for t ∈ [0, 1).
With this strategy, the existence of zeros of ǫ∞ r,n is no longer an
issue; for any r ≥ 0, simply examine the finite series approximation
of |ǫr,n (z(t))|, 0 ≤ t < 1. It turned out in practice, however, that for
each n examined, setting r to be the largest zero of ǫ∞ r,n produced the
smallest uniform error bound for |ǫr,n (z(t))|.

8.5.2 Zeros of ǫ∞
r,n

The question remains, for n fixed, does ǫ∞ r,n have zeros as a function of
r? The answer was found to be yes in all cases examined. Refer to Fig-
ure 8.1 for plots of [ǫ∞
r,n ]
1/5
, n = 0, 1, 2, 3, showing the location of zeros
2
for the first few cases. The plots of Figure 8.1 motivated an extensive
numerical investigation resulting in empirical data on the number, the
smallest and the largest zeros of ǫ∞ r,n n = 0, . . . , 60. The results of this
investigation are listed in Tables C.1 and C.2 of Appendix C. From
these tables, ǫ∞r,n appears to have about 2n zeros, and the largest zero
appears to be about size n. Beyond these superficial observations, there
is no more obvious trend in the data.

8.5.3 Improved Lanczos Formula

As an example of how these observations can be applied to achieve
improved error bounds, consider the case n = 6 which is the last case
for which Lanczos gives an explicit formula and error bounds in [14].
2 The
functions ǫ∞
r,n decrease rapidly with n, yet grow very quickly once r ≈ n,
and so to show the plots on a single set of axes the fifth roots of these
functions are plotted.

111
Chapter 8. Implementation & Error Discussion

1/5
[ǫ∞
r,1 ]

1/5
[ǫ∞
r,3 ]

r
5

1/5
[ǫ∞
r,2 ]

1/5
[ǫ∞
r,0 ]

−1

1/5
Figure 8.1: ǫ∞
r,n , n = 0, 1, 2, 3, −1/2 < r < 6

That is,
6
" #
√ X ′
Γ(z+1) = 2π (z+r+1/2)z+1/2 e−(z+r+1/2) ak (r)Hk (z) + ǫr,6 (z) .
k=0

With the choice of r = 5, Lanczos gives a uniform error bound of

|ǫr,6 (z)| < 2 × 10−10 .
Twelve zeros of ǫ∞
r,6 were found; call these r0 , r1 , . . . , r11 . For each
of these zeros, using
the estimation techniques of Section 8.5.1, the
maximum Mrj ,6 of ǫrj ,6 (it/(1 − t)) was estimated along with its loca-
tion on the imaginary axis using first fifteen terms of ǫrj ,6 (it/(1 − t)).

112
Chapter 8. Implementation & Error Discussion

.
Table 8.4 summarizes the results. Notice that the largest zero, r =
6.779506 yields the least maximum error of 2.72 × 10−12 , nearly 1/100
of Lanczos’ estimate.

j rj t/(1 − t) Mrj ,6
0 −0.117620 0.565599 4.71 × 10−4
1 0.684391 1.525502 2.75 × 10−6
2 1.450013 2.388058 8.88 × 10−8
3 2.182290 3.172714 6.78 × 10−9
4 2.883225 3.887447 9.30 × 10−10
5 3.553321 4.538907 1.99 × 10−10
6 4.191832 5.133674 6.07 × 10−11
7 4.796781 5.678871 2.49 × 10−11
8 5.364813 6.183352 1.30 × 10−11
9 5.891184 6.660957 8.02 × 10−12
10 6.372580 7.137483 5.29 × 10−12
11 6.779506 7.883760 2.72 × 10−12

Table 8.4: Approximate Zeros of ǫ∞

r,6

Plots of log ǫrj ,6 (it/(1 − t)) for j = 0, 1, . . . , 11 are shown in Fig-
ure 8.2. Observe that the maximum of the curve associated with the
largest zero r11 gives the smallest maximum error of approximately
exp (−26.6).
To get an idea of how the maximum error along the imaginary axis
varies with r , refer to Figure 8.3 in which log |ǫr,4 (it/(1 − t))| is plotted
for 0 < t < 1, 0 < r < 5. The ends of the deep troughs near t = 1
correspond to the zeros of ǫ∞ 3
r,4 . The projection of the troughs onto
the t − r plane are not lines parallel to the t-axis, but rather slowly
varying curves. Bearing in mind that the vertical scale on the plot is
logarithmic, the effect of the the r parameter on the error is dramatic.
3 Forgraphing purposes the t range was truncated slightly before t = 1, and
so the vertical asymptotes associated with the ends of the troughs are not
present in the plot.

113
Chapter 8. Implementation & Error Discussion

0 t 1

−7 r0

r1

r2
r3
−21

r11

−35

Figure 8.2: log ǫrj ,6 (it/(1 − t)), 0 < t < 1, j = 0, . . . , 11

8.5.4 Optimal Formulas

In practice, for each value of n examined, the largest zero of ǫ∞ r,n was
found to produce the smallest maximum error (among all zeros), and
bounds much improved on Lanczos’ estimates. Tables C.1 and C.2 of
Appendix C lists, for n = 0, . . . , 60, the largest zero of ǫ∞ r,n , denoted
r(n), the corresponding estimates
Mr(n),n,5 and Mr(n),n,15 of the maxi-
mum of ǫr(n),n (it/(1 − t)) using five and then fifteen terms

of the sum,
respectively, and finally, the location of the maximum of ǫr(n),n (it) on
the imaginary
axis.4 . From
this data it appears that five terms of the
series ǫr(n),n (it/(1 − t)) suffice to estimate the maximum error in all
4 Ofparticular interest is the n = 1 case. In Lanczos’ original paper [14] he
remarks:

Particularly remarkable is the approximation of only two terms

114
Chapter 8. Implementation & Error Discussion

−6

−22

0
t

1 r
5

Figure 8.3: log |ǫr,4 (it/(1 − t))|, 0 < t < 1, 0 < r < 5

but the first few cases.

Observe that the n = 10, r = 10.900511 row of Table C.1 is sufficient
to guarantee 16 digit floating point accuracy in the right-half plane. In
this case, resolving the series Sr(10),10 into partial fractions and rescaling
(γ = 1.5):
√

z+1/2 −(z+2) 1.084635
z! = (z + 2) e 2π 0.999779 + ,
z+1

which is correct up to a relative error of 2.4 · 10−4 everywhere

in the right complex half plane.

Note that we use r to denote what Lanczos called γ to avoid confusion with
Euler’s constant. It turns out that for n = 1, r = 1.5 is coincidentally very
close to 1.489193 · · · , the largest zero of ǫ∞
r,1 .

115
Chapter 8. Implementation & Error Discussion

the coefficients as in equation (6.14) yields the approximation

z+1/2 " 10
r #
e z + r(10) + 1/2 X dk
Γ(z + 1) ≈ 2 d0 + ,
π e k=1
z+k

where the coefficients dk are as given in Table 8.5.

k dk
0 +2.48574089138753565546 × 10−5
1 +1.05142378581721974210 × 100
2 −3.45687097222016235469 × 100
3 +4.51227709466894823700 × 100
4 −2.98285225323576655721 × 100
5 +1.05639711577126713077 × 100
6 −1.95428773191645869583 × 10−1
7 +1.70970543404441224307 × 10−2
8 −5.71926117404305781283 × 10−4
9 +4.63399473359905636708 × 10−6
10 −2.71994908488607703910 × 10−9

Table 8.5: Sample coefficients, n = 10, r = 10.900511

8.6 Additional Remarks

Although the results of Section 8.5 are empirical in nature, the resulting
formulas and corresponding error bounds seem promising, and lead to
many unanswered questions. In particular,

(i) For each n ≥ 0, how many zeros does ǫ∞

r,n have? What is the
largest zero?
(ii) For fixed n, why does the largest zero of ǫ∞
r,n seem to give rise to
the smallest uniform bound on |ǫr,n (z)|?
(iii) For fixed n and r, is it possible to analytically determine the
maximum of |ǫr,n (it/(1 − t))| for t ∈ [0, 1)?
(iv) Perhaps most importantly, given ǫ > 0, what is the best choice of
r and n to estimate Γ with a relative error of at most ǫ?

116
Chapter 8. Implementation & Error Discussion

It was observed early in the investigation that, not only does ǫ∞r,n
have zeros as a function of r, but so too do the individual ak (r). This
is not surprising, since
∞
X
ǫ∞
r,n = ak (r) ≈ an+1 (r) ,
k=n+1

assuming the coefficients are rapidly decreasing. However, the relation-

ship between ǫ∞ r,n and an+1 (r) is eye opening; refer to Figures 8.4 and
8.5 for plots of [ǫ∞
r,6 ]
1/7
and [a7 (r)]1/7 , respectively. Both functions are
then plotted in Figure 8.6 and the plots are virtually indistinguishable.
0.5

[ǫ∞
r,6 ]
1/7

5
r

−0.5

Figure 8.4: [ǫ∞

r,6 ]
1/7
, −0.2 < r < 7

Numerical checks show that the zeros of ǫ∞ r,n and an+1 (r) nearly co-
incide in the cases examined. Furthermore, at r(n) equal the largest
zero of ǫ∞
r,n , an+1 (r(n)) appears to be a very good estimate for the max-
imum of |ǫr,n (it/(1 − t))| on [0, 1). This observation seems in-keeping
with the rule of thumb for asymptotic series with rapidly decreasing
terms: the error is about the size of the first term omitted. Oddly,
though, this is not the situation here: the second coefficient omitted is
almost exactly the same size as the first, but of opposite sign. These
observations are summarized in Table 8.6 where for each n = 0, . . . , 12,

117
Chapter 8. Implementation & Error Discussion
0.5

[a7 (r)]1/7

5
r

−0.5

Figure 8.5: [a7 (r)]1/7 , −0.2 < r < 7

0.5

5
r

−0.5

Figure 8.6: [ǫ∞

r,6 ]
1/7
and [a7 (r)]1/7 , −0.2 < r < 7

listed are r(n), the estimated uniform error bound Mr(n),n,15 , and the
first two terms omitted from the series, an+1 (r(n)) and an+2 (r(n)).

118
Chapter 8. Implementation & Error Discussion

n r(n) Mr(n),n,15 an+1 (r(n)) an+2 (r(n))

0 0.319264 5.5 × 10−3 +5.4 × 10−3 −7.7 × 10−3
1 1.489194 1.0 × 10−4 −1.0 × 10−4 +1.1 × 10−4
2 2.603209 6.3 × 10−7 +5.3 × 10−7 −3.4 × 10−7
3 3.655180 8.5 × 10−8 +8.4 × 10−8 −9.3 × 10−8
4 4.340882 4.3 × 10−9 +4.2 × 10−9 −4.6 × 10−9
5 5.581000 1.2 × 10−10 −1.2 × 10−10 +1.2 × 10−10
6 6.779506 2.7 × 10−12 +2.7 × 10−12 −2.5 × 10−12
7 7.879012 3.9 × 10−14 +3.6 × 10−14 −4.7 × 10−14
8 8.406094 6.9 × 10−15 +6.9 × 10−15 −7.1 × 10−15
9 9.656578 2.1 × 10−16 −2.0 × 10−16 +2.0 × 10−16
10 10.900511 6.1 × 10−18 +6.1 × 10−18 −5.9 × 10−18
11 12.066012 1.1 × 10−19 −1.1 × 10−19 +9.1 × 10−20
12 13.144565 5.2 × 10−21 −5.1 × 10−21 +5.6 × 10−21

Table 8.6: Comparison of Mr(n),n and an+1 (r(n))

The data in Table 8.6 would lead one to conjecture that

Mr(n),n ≈ |an+1 (r(n))| ,

which, if true, provides a much simpler method for bounding |ǫr(n),n (z)|
in the right-half plane.

119
Chapter 9
Comparison of the Methods

In this chapter we revisit the computation of Γ(1 + z) from Chapter 2

for the purpose of comparing Lanczos’ method against those of Stirling
and Spouge. For each of these methods Γ(20 + 17i) is computed with a
relative error |ǫρ | < 10−32 , and then formulas for Γ(1 + z) are given for
each with a uniform error bound of 10−32 in the entire right-half plane
Re(z) ≥ 0.

9.1 Stirling’s Series

For this section recall the notation of Section 2.5: let
Z ∞
B2n (x)
EN,n (z) = dx ,
0 2n(z + N + x)2n
denote the integral error term of equation (2.11), and denote by UN,n (z)
the upper bound on EN,n (z) given by Theorem 2.3:
2n+2
1 B2n+2
UN,n (z) = .
cos (θ/2) (2n + 2)(2n + 1)(z + N)2n+1

9.1.1 Γ(20 + 17i) with |ǫρ | < 10−32

Since z = 19 + 17i is fixed, for convenience write EN,n = EN,n (19 + 17i)
and UN,n = UN,n (19 + 17i).
From equation (8.2), EN,n is related to the relative error ǫρ via
1
eEN,n = ,
1 − ǫρ

120
Chapter 9. Comparison of the Methods

from which
ǫρ = 1 − e−EN,n .
We want |ǫρ | < 10−32 which means |1 − e−EN,n | < 10−32 , which is
guaranteed if

|EN,n | ≤ |UN,n | < 10−32 /e ≈ 3.6 × 10−33 .

The pairs (N, n) which achieve this bound are listed in Table 9.1.

N n |UN,n | <
18 10 3.5 × 1033
10 11 2.9 × 1033
4 12 3.3 × 1033
0 13 2.6 × 1033

Table 9.1: Upper bound on |EN,n (19 + 17i)| in Stirling Series

Selecting (N, n) = (4, 12), equation (2.11) yields

" 4
#
Y
log Γ(20 + 17i) (z + k)
k=1
1
≈ (19 + 17i + 4 + 1/2) log (19 + 17i + 4) − (19 + 17i + 4) + log 2π
2
12
X B2j
+ ,
j=1
2j(2j − 1)(19 + 17i + 4)2j−1

from which
.
Γ(20 + 17i) = −6.6530978807100357093202320786706 × 1013
+ i 1.3813486137818296429873066956513 × 1014 .

Comparing this value against a high precision value produced by Maple,

we find a relative error of approximately 8.5 × 10−34 .
.
The Bernoulli numbers in the sum range in size from B2 = 1.7×10−1
.
to B24 = 8.7 × 104 . The corresponding terms of the sum range in size
from 2.9 × 10−3 to 5.0 × 10−32 .

121
Chapter 9. Comparison of the Methods

9.1.2 Γ(1 + z) with |ǫρ | < 10−32 Uniformly

The goal this time is a formula for Γ(1 + z) with relative error |ǫρ |
bounded uniformly by 10−32 in the right-half plane. From the analysis
in Section 2.5.2, |UN,n (z)| is maximized at z = 0, and N must be at
least 1 for otherwise |UN,n (z)| becomes unbounded as z → 0+ along
the real axis. By the same argument as in the previous section, again
in this case the bound |UN,n (0)| < 3.6 × 10−33 is sufficient to ensure
|ǫρ | < 10−32 . Computing |UN,n (0)| for pairs (N, n) which achieve this
bound results in Table 9.2. Selecting (N, n) = (17, 17), for example, we

N n |UN,n (0)| <

19 15 3.5 × 10−33
18 16 1.4 × 10−33
17 17 9.4 × 10−34
16 18 9.7 × 10−34
15 19 1.7 × 10−33
14 21 1.1 × 10−33
13 23 2.4 × 10−33
12 28 2.2 × 10−33

Table 9.2: Upper bound on |EN,n (0)| in Stirling Series

thus obtain the following formula

" 17
#
Y
log Γ(1 + z) (z + k)
k=1
1
≈ (z + 17 + 1/2) log (z + 17) − (z + 17) + log 2π
2
17
X B2j
+ .
j=1
2j(2j − 1)(z + 17)2j−1

In this case, the Bernoulli numbers in the sum range in size from
. .
B2 = 1.7 × 10−1 to B34 = 4.3 × 1011 , while the corresponding terms of
the sum with z = 0 range in size from 4.9 × 10−3 to 9.5 × 10−33 .

122
Chapter 9. Comparison of the Methods

9.2 Spouge’s Method

Before we begin, note
√ the simplification of equation (2.17) given by the
cancellation of the 2π exp(−a) term which occurs in the leading fac-
tor, and (2π)−1/2 exp(a) which occurs in each of the coefficients (2.18).
Spouge’s formula may then be stated
" N
#
√ X d k (a)
Γ(z + 1) = (z + a)z+1/2 e−z 2π e−a + + ǫ(z) (9.1)
z+k
k=1

where N = ⌈a⌉ − 1, and

(−1)k−1
dk (a) = (−k + a)k−1/2 e−k . (9.2)
(k − 1)!
Under this simplification, the relative error in Spouge’s formula remains
the same: √
a 1
|ǫS (a, z)| < .
(2π)a+1/2 Re(z + a)

9.2.1 Γ(20 + 17i) with |ǫρ | < 10−32

For |ǫρ | < 10−32 we require
√
a 1
a+1/2
< 10−32 .
(2π) Re(19 + 17i + a)
The left hand side of this expression is eventually decreasing as a func-
tion of a, and drops below 10−32 between a = 38 and a = 39. Thus
N = 38, and computing the coefficients using (9.2) results in Table 9.3.
Inserting these values in equation (9.1) we find to 32 digits
.
Γ(20 + 17i) = −6.6530978807100357093202320786706 × 1013
+ i 1.3813486137818296429873066956513 × 1014 .

In this instance, comparing this value against a high precision value

produced by Maple, the relative error is approximately 2.6 × 10−44 ,
which is much more than necessary. Through experimentation the value
a = 29 is found give the desired accuracy, and the series can thus be
reduced by 10 terms.

123
Chapter 9. Comparison of the Methods

k dk (39) k dk (39)
0 +2.8947105233858572256681259383826×10−17 20 −1.4612043596445044389340540663592×10−1
1 +2.2677611785616408448051623649413×100 21 +1.6856896402284750130501640213058×10−2
2 −3.0458858426261684740752182066857×101 22 −1.5553541637142526870805344933539×10−3
3 +1.9357212181425501030368331204744×102 23 +1.1302120506210549847428847358464×10−4
4 −7.7429949716371535150533416191734×102 24 −6.3472060001627477262571023780708×10−6
5 +2.1876179361012802304828398825556×103 25 +2.6919634338361067149994936318051×10−7
6 −4.6438537277548560417623360405775×103 26 −8.3801786194339100406682820681103×10−9
7 +7.6927388480783573540015702664876×103 27 +1.8481319798192137441154507903452×10−10
8 −1.0195917484892786070501293142239×104 28 −2.7610263895691596705265848094462×10−12
9 +1.0999167217153749462567889811423×104 29 +2.6382697777056113806964677703858×10−14
10 −9.7740178567227019446568421903035×103 30 −1.4954804884394724082929843529081×10−16
11 +7.2136820304845623606047714747868×103 31 +4.5441806633655985113201441255734×10−19
12 −4.4462490659018128337839056551451×103 32 −6.4289947451722952263534424159053×10−22
13 +2.2961572174783000519663098692718×103 33 +3.4516439287785093504202451110646×10−25
14 −9.9491746363493229520431441706026×102 34 −5.1380351805226988204329822384216×10−29
15 +3.6160749560761185688455782354331×102 35 +1.2606607461787976943234811185446×10−33
16 −1.1004495488313215123003718194264×102 36 −1.9452281496812994762288072286241×10−39
17 +2.7947852816195578994204731195154×101 37 +2.2292761113182594463394149297371×10−47
18 −5.8947917688963354031891858182660×100 38 −2.2807244297698558308437038776471×10−60
19 +1.0259323949497584239156793246602×100

Table 9.3: Coefficients of Spouge’s series, a = 39

9.2.2 Γ(1 + z) with |ǫρ | < 10−32 Uniformly

In the case of Spouge’s method with a uniform bound the procedure is
essentially the same. The bound on the relative error is greatest when
Re(z) = 0, so for |ǫρ | < 10−32 we require
√
a 1
a+1/2
< 10−32 .
(2π) a

Again this is satisfied for a between 38 and 39, so N = 38 and the

coefficients in the uniform case are those of Table 9.3.

124
Chapter 9. Comparison of the Methods

9.3 Lanczos’ Method

In the Lanczos case there is only one situation to consider since the error
bounds developed in Chapter 8 give only uniform error bounds. For
this method we examine pairs (n, r(n)), where r(n) is the largest zero of
ǫ∞
r,n , and determine the least n which yields the prescribed error bound.
For |ǫρ | < 10−32 we find using Table C.1 that n = 21, r = 22.618910
gives the empirical bound |ǫr,n (it/(1 − t))| ≤ 2 × 10−34 for t ∈ [0, 1).
Thus, by equation (8.6),
π 1/2
|ǫρ | ≤ 2 × 10−34
e

< 2.2 × 10−34 .

Using the formulation (6.14), we find

r z+1/2 " 21
#
e z + r(21) + 1/2 X dk
Γ(z + 1) = 2 d0 +
π e k=1
z+k

where the coefficients dk are given in Table 9.4.

Evaluating this expression at z = 19 + 17i we find to 32 digits
.
Γ(20 + 17i) = −6.6530978807100357093202320786706 × 1013
+ i 1.3813486137818296429873066956513 × 1014 .

The relative error in this calculation is less than 1.5 × 10−42 .

9.4 Discussion
All three methods considered here have their benefits and shortfalls,
and the question of which is best is not a clearcut one. Several factors
must be considered, in particular:

1. The computational cost of the method;

2. Whether the series coefficients are available in precomputed form;

125
Chapter 9. Comparison of the Methods

k dk
0 +2.0240434640140357514731512432760 × 10−10
1 +1.5333183020199267370932516012553 × 100
2 −1.1640274608858812982567477805332 × 101
3 +4.0053698000222503376927701573076 × 101
4 −8.2667863469173479039227422723581 × 101
5 +1.1414465885256804336106748692495 × 102
6 −1.1135645608449754488425056563075 × 102
7 +7.9037451549298877731413453151252 × 101
8 −4.1415428804507353801947558814560 × 101
9 +1.6094742170165161102085734210327 × 101
10 −4.6223809979028638614212851576524 × 100
11 +9.7030884294357827423006360746167 × 10−1
12 −1.4607332380456449418243363858893 × 10−1
13 +1.5330325530769204955496334450658 × 10−2
14 −1.0773862404547660506042948153734 × 10−3
15 +4.7911128916072940196391032755132 × 10−5
16 −1.2437781042887028450811158692678 × 10−6
17 +1.6751019107496606112103160490729 × 10−8
18 −9.7674656970897286097939311684868 × 10−11
19 +1.8326577220560509759575892664132 × 10−13
20 −6.4508377189118502115673823719605 × 10−17
21 +1.3382662604773700632782310392171 × 10−21

Table 9.4: Coefficients of Lanczos’ series, r = 22.618910

3. The precision of the available hardware and software; and

4. Whether a uniform or point-wise error bound is required.
All three methods considered here have essentially the same leading
factor but differ in the terms of their respective series. Assuming for the
moment that the constant terms of the series have been precomputed,
the computational cost of the series terms then serves as a measure
of the efficiency of each method. Stirling’s series has the additional
product term which must also be considered in the cost.
For the evaluation of Γ(1 + z) with a uniformly bounded error, the
Lanczos method appears to be the clear winner. In the case |ǫρ | < 10−32
considered here, 21 divisions and 21 additions were required in the

126
Chapter 9. Comparison of the Methods

sum. By contrast, Stirling’s series required 16 powers, divisions and

additions in the series, and a further 16 multiplications in the product
term. Spouge’s method requires 38 additions and divisions.
For evaluation at a single value of z, especially if |z| is large, Stir-
ling’s series is most efficient due to the error term which decreases like
|z+N|−(2n+1) with increasing |z|. Spouge’s error decreases with increas-
ing |z|, but only like |Re(z + a)|−1 . The Lanczos error also decrease to
zero as |z| → ∞, by virtue of the choice of r as the largest zero of ǫ∞ r,n ,
but the rate has not been quantified, though experimentally it appears
to be slow.
The situation changes, however, if the constant terms in the series
must be computed first. In this case, the easiest coefficients to evaluate
are those of Spouge’s series as given by equation (9.2). Simple recursion
formulas for the calculation of the Bernoulli numbers in Stirling’s series
are known, see [28, p.6], which puts Stirling’s series in second place. The
coefficients of Lanczos’ series are the most labour intensive to calculate.
In the examples presented here the focus was relative error, and
we have tacitly assumed infinite floating point precision, an unrealis-
tic assumption in practice. In reality, the floating point precision of
the hardware and software used may introduce round-off or overflow
errors, and this becomes a serious concern in gamma function calcu-
lations. To prevent overflow error in Spouge’s and Lanczos’ methods,
one should first compute log Γ(1 + z) as in Stirling’s series, and then
take the exponential of the result.
As was already noted in Chapter 6, generous decimal precision is
required in the calculation of coefficients in the Lanczos method since
the relatively small coefficients are the result of adding many larger
terms of alternating sign. The matrix methods of Section 6.7 help
to lessen this problem. In Stirling’s series one must pay attention to
the rapidly increasing Bernoulli numbers and select the parameters n
and N accordingly. As was noted in Section 2.5.2, it is possible to
achieve a relative error of order 10−34 in the calculation of Γ(7 + 13i)
using (N, n) = (0, 38), but the largest Bernoulli number required will
be of order 1050 . Roundoff can also be a problem in the summation of
the series themselves, even with accurately computed coefficients. For
example, the coefficients in Table 9.3 span 63 orders of magnitude, and
many significant digits of some individual terms will be lost in a 32
digit environment.

127
Chapter 10
Consequences of Lanczos’
Paper

There are several consequences to the techniques and ideas used in

Lanczos’ paper which extend beyond the computation of the gamma
function. This chapter illustrates some of these extensions and men-
tions areas worthy of further study.
We first see how Lanczos’ derivation of (1.1) is a special case of a
much more general process which defines a transform between square
summable functions on [−π/2, π/2] and analytic functions on half planes.
This can be used to express certain functions as series similar to (1.1).
Some less than obvious combinatorial identities are then stated as con-
sequences of the Lanczos Limit Formula. Finally, an improvement to
Stirling’s formula is noted based on the existence of zeros of the error
function ǫ∞r,0 .

10.1 The Lanczos Transform

In his original paper, Lanczos refers to equation (3.6) as an integral
transform, which is here termed a “Lanczos Transform”. √ That is,
−(z+1/2)
in (3.6), Γ(z +1/2)(z +r +1/2) exp (z + r + 1/2) 2 is the Lanc-
zos Transform of fE,r (θ). The transform is worthy of special attention
as the techniques used to compute it generalize easily from those used
in the gamma function case.
Suppose F (z) is defined for Re(z) ≥ 0 as
Z π/2
F (z) = cos2z θg(θ) dθ
−π/2

128
Chapter 10. Consequences of Lanczos’ Paper

where g(θ) ∈ L2 [−π/2, π/2] is even. Note that g is very general, the
only requirement being that it be square summable with respect to
Lebesgue measure. Then the Fourier coefficients

2 π/2
Z
ak = g(θ) cos (2kθ) dθ .
π −π/2

are defined, and we can associate with g its Fourier series

∞
X ′
g(θ) ∼ ak cos (2kθ) .
k=0

The ∼ symbol is used since the series converges to g in L2 [−π/2, π/2],

though not necessarily pointwise. Nonetheless, since cos2z θ ∈ L2 [−π/2, π/2],
it follows from Parseval’s theorem and identity (3.9) of Lemma 3.2 that
Z π/2
F (z) = cos2z θg(θ) dθ
−π/2

∞
√ Γ(z + 1/2) X′
= π ak Hk (z) . (10.1)
Γ(z + 1) k=0

Each of the functions Hk (z) is analytic for Re(z) > −1/4. Let K
be a compact subset of this region. Taking r = 1/4 in the proof of
Theorem 5.1 shows that Hk (z) = O k −1/2−δ on K for some δ > 0. By
the Schwarz inequality,
∞
!2 ∞
! ∞ !
X X X 1
|ak ||Hk (z)| ≤ C |ak |2
k=N k=N k=N
k 1+2δ

on K for some constant C which does not depend on N. The right hand
side of this inequality can be made as small as desired by choosing N
sufficiently large. Thus the series in (10.1) converges absolutely and
uniformly on K, and since K was arbitrary, (10.1) defines F (z) as an
analytic function on Re(z) > −1/4.
If in addition g (k) (θ) is continuous with g (k) (−π/2) = g (k) (π/2),
k = 0, . . . , m, m ≥ 1, then the Fourier coefficients ak = O (k −m ). For

129
Chapter 10. Consequences of Lanczos’ Paper

an arbitrary compact K ⊂ C \ {−1, −2, . . .}, again from the proof of

Theorem 5.1 we have Hk (z) = O (k −2σ−1 ) on K. Therefore, provided
Re(−2z − m − 1) < −1, the series in (10.1) will converge absolutely
and uniformly on K, and thus defines an analytic function in the region

{z ∈ C | Re(z) > −m/2 and z 6= −1, −2, −3, . . .} .

Since the Γ(z+1/2)/Γ(z+1) factor of (10.1) cancels the poles at the neg-
ative integers but introduces simple poles at z = −1/2, −3/2, −5/2, . . .,
the conclusion is that equation (10.1) defines F (z) as an analytic func-
tion on the half plane Re(z) > −m/2 except for z = −1/2, −3/2, −5/2, . . ..
Furthermore, if F (z) is approximated by truncating
P∞ the series in (10.1)
th
after the n term, the resulting error ǫr,n (z) = k=n+1 ak Hk (z) has the
property that in the half plane Ω = Re(z) ≥ 0,
∞
X
lim ǫr,n (z) = ak
|z|→∞
z∈Ω k=n+1

n
X
= g(0) − ak ,
k=0

a constant. Invoking Lemma 8.2, the maximum of |ǫr,n (z)| in Ω occurs

on the line Re(z) = 0, and the maxima can be located and estimated
by examining the first few terms of series |ǫr,n (it/(1 − t))|, 0 ≤ t < 1.
As in the gamma function case, the coefficients can be found in
several ways. If F (z) is known on the non negative integers, the
Horner type method of Section 6.2, and the Lanczos method of Sec-
tion 6.1 carry over word for word, with the Horner method probably
the more practical. If g(θ) is easy to compute, the finite Fourier se-
ries described in Section 6.4.1 may prove more efficient. The Lanczos
method is repeated here. From the definition of Chebyshev polynomi-

130
Chapter 10. Consequences of Lanczos’ Paper

als, Tk (cos θ) = cos (kθ). Thus

π/2
2
Z
ak = g(θ) cos (2kθ) dθ
π −π/2

π/2 k
2
Z X
= g(θ) C2j,2k cos2j θ dθ
π −π/2 j=0

k π/2
2X
Z
= C2j,2k cos2j θg(θ) dθ
π j=0 −π/2

k
2X
= C2j,2k F (j) .
π j=0

To go in the other direction, suppose the function F (z) is given

analytic on a domain including the non-negative real axis. Compute
the coefficients
k
2X
ak = C2j,2k F (j)
π j=0

and form the formal sum

∞
X ′
g(θ) = ak cos (2kθ) .
k=0
P∞
If this series is in L2 [−π/2, π/2] , that is, if 2
k=0 ak < ∞, then it
defines g(θ) as the inverse transform of F (z).

An Example
The following example is found in [13, p.45]. The Bessel function of
the first kind of order z not necessarily an integer can be expressed
Z π/2
1 (x/2)z
Jz (x) = √ cos2z θ cos (x sin θ) dθ .
π Γ(z + 1/2) −π/2

131
Chapter 10. Consequences of Lanczos’ Paper

Letting √
F (z, x) = π Jz (x)Γ(z + 1/2)(x/2)−z ,
we see that F (z, x) is the Lanczos transform of cos (x sin θ). Hence
∞
√ Γ(z + 1/2) X′
F (z, x) = π ak (x)Hk (z)
Γ(z + 1) k=0

so that ∞
(x/2)z X′
Jz (x) = ak (x)Hk (z) .
Γ(z + 1) k=0
The coefficients ak (x) are given by
k
2X
ak (x) = C2j,2k F (j, x)
π j=0

k
2 X
=√ C2j,2k Jj (x)Γ(j + 1/2)(x/2)−j .
π j=0

On the other hand, with the transformation φ = θ−π/2, the Fourier

coefficients become

2 π/2
Z
ak (x) = cos (x sin θ) cos (2kθ) dθ
π −π/2

π
2
Z
= cos (x sin(φ + π/2)) cos (2k(φ + π/2)) dφ
π 0

π
2
Z
k
= (−1) cos (x cos φ) cos (2kφ) dφ
π 0

= 2J2k (x) .

This shows how the Bessel function of even integer order 2k can be
expressed in terms of Bessel functions of integer order less than or

132
Chapter 10. Consequences of Lanczos’ Paper

equal to k:
k
1 X
J2k (x) = √ C2j,2k Jj (x)Γ(j + 1/2)(x/2)−j .
π j=0

More importantly,
∞
2(x/2)z X
Jz (x) = J2k (x)Hk (z) .
Γ(z + 1) k=0

It is interesting in this case to examine what the properties of g(θ) =

cos (x sin θ) say about the rate of convergence of the series. In this
case, g (k) (−π/2) = g (k) (π/2), for all k ≥ 0, so that for fixed x, ak (x) =
2J2k (x) = O (k −n ) for any n ≥ 0 as k → ∞. Furthermore, again for
fixed x, Jz (x) is entire as a function of z.

10.2 A Combinatorial Identity

The following non-trivial identities are consequences of the Lanczos
Limit Formula 7.1:
Lemma 10.1. Suppose n ≥ 1 is an integer. Then
(i)
n
1 X n!n!
+ (−1)k =0
2 k=1 (n − k)!(n + k)!

(ii)
n
X k 2n 1
(−1)n+k =
k=1
(n − k)!(n + k)! 2

Proof of Lemma 10.1: From 7.1,

" ∞
#
z 1 k −k 2 /r
X
Γ(z + 1) = 2 lim r + (−1) e Hk (z) .
r→∞ 2
k=1

Setting z = n terminates the series giving

n
" #
1 X 2 n!n!
n! = 2 lim r n + (−1)k e−k /r .
r→∞ 2 (n − k)!(n + k)!
k=1

133
Chapter 10. Consequences of Lanczos’ Paper

Since r n → ∞ while n! remains constant, it follows that the term in

square brackets must tend to zero as r → ∞. This gives (i).
For (ii), again setting z = n terminates the series which may be
written
+ nk=1 (−1)k e−k /r (n−k)!(n+k)!
1
P 2 n!n!
2
n! = 2 lim .
r→∞ r −n
Applying L’Hôpital’s rule n times then gives
Pn k −k 2 /r n!n!
H k=1 (−1) e (n−k)!(n+k)!
k2
n! = 2 lim
r→∞ −nr −n+1
Pn k −k 2 /r n!n!
H k=1 (−1) e (n−k)!(n+k)!
k4
= 2 lim
r→∞ n(n − 1)r −n+2
..
.
Pn k −k 2 /r n!n!
H k=1 (−1) e (n−k)!(n+k)!
k 2n
= 2 lim
r→∞ (−1)n n!
n
X n!
= 2 (−1)k+n k 2n
k=1
(n − k)!(n + k)!

which yields (ii) upon canceling the n! terms and moving the constant
2 to the left hand side.

10.3 Variations on Stirling’s Formula

For large real z > 0, Stirling’s formula (2.12) is often expressed
Γ(z + 1) = e−z z z+1/2 (2π)1/2 eθ/(12z) (10.2)
for some 0 < θ < 1, with resulting relative error exp (θ/(12z)) − 1 ≈
θ/(12z) which decreases with increasing z. However, without some
manipulation, the formula is not very effective for small z. Spouge
in [27] derives the more visually pleasing approximation
√
Γ(z + 1) ≈ (z + 1/2)z+1/2e−(z+1/2) 2π (10.3)
with bounded relative error
1 log 2 1
ǫ(z) < √ ,
2πe π σ + 1/2

134
Chapter 10. Consequences of Lanczos’ Paper

which is at most 0.107 at σ = Re(z) = 0. As he notes, (10.3) is not

only simpler, but more accurate than Stirling’s formula (10.2).
In this section a formula similar in form to both (10.2) and (10.3)
is found as a result of the previous analysis of the relative error in
Section 8.5.4, and a generalization of the formula is noted. For this
section the variable a := r + 1/2 is introduced to simplify the resulting
formulas.
The work in Section 8.5.4 furnishes an improvement to (10.2), albeit
with the introduction of a slightly complicated constant. Letting r(0)
denote the largest zero of ǫ∞
r,0 , and setting a = r(0) + 1/2, we find

√

(z+1/2) −(z+a) a0 (r(0))
Γ(z + 1) = (z + a) e 2π + ǫr(0),0 (z)
2

√
≈ (z + a)(z+1/2) e−(z+a) 2π (10.4)

since ǫ∞ r(0),0 = 0 implies a0 (r(0))/2 = 1. This approximation is valid

in the entire right half plane Re(z) ≥ 0, is asymptotic to Γ(z + 1) as
|z| → ∞ in this region, and has an empirically determined relative error
|ǫr(0),0 (z)| < 0.006.
.
The constant a = 0.819264 appears complicated, but it turns out
that the solutions of ǫ∞ r,n = 0 can be found explicitly in terms of Lambert
W functions in the n = 0 case, as will now be shown.

Lemma 10.2. Let r(0) be the largest zero of ǫ∞

r,0 , and a = r(0) + 1/2.
Then a = −W (−1/π)/2.

Proof of Lemma 10.2: Recall that

a0 er+1/2
=p
2 2π(r + 1/2)

ea
=√ .
2πa

135
Chapter 10. Consequences of Lanczos’ Paper

The real zeros of ǫ∞

r,0 are the solutions of 1 − a0 /2 = 0, so

ea
1= √
2πa

2πa = e2a

1
−2ae−2a = − .
π

That is, −2a = W (−1/π) whose only real roots are

a = −W0 (−1/π)/2

.
= 0.276914

and
a = −W−1 (−1/π)/2

.
= 0.819264 .

The value of a = −W−1 (−1/π)/2 makes (10.4) exact at z = 0,

and (10.4) is exact at infinity, regardless of the value of a. It is natural
to ask if, given z, can a value for a be found such that (10.4) is again
exact? In other words, can a function a(z) be defined on some subset
of C such that on this set
√
Γ(z + 1) = (z + a(z))(z+1/2) e−(z+a(z)) 2π ?
Recalling the definition of Fr (z) from Section 8.2
Fr (z) = Γ(z + 1)(z + r + 1/2)−(z+1/2) ez+r+1/2 (2π)−1/2 ,
for fixed z ≥ 0 real, this is equivalent to asking if Fr (z) − 1 = 0 has
solutions. The answer is yes. The following property of Fr (z) will be
needed (see [27]):

136
Chapter 10. Consequences of Lanczos’ Paper

Lemma
p 10.3. F0 (z) strictly increases to 1 as z → ∞, and F0 (0) =
e/π .

From this result follows

Lemma 10.4. Let z > 0 be fixed. Then Fr (z) − 1 = 0 has solutions as
a function of r .

Proof of Lemma 10.4: First, note that Fr (z) is continuous for (r, z) ∈
(−1/2, ∞) × [0, ∞). Secondly, observe that

lim Fr (z) − 1 = lim Fr (z) − 1 = ∞ .

r→−1/2 r→∞

Thirdly, by Lemma 10.3, F0 (z) −1 ↑ 0 as z → ∞, so that F0 (z) −1 < 0.

Thus for z > 0 fixed, Fr (z) has at least two zeros as a function of r.

An explicit formula for a(z) can be found with the help of the Lam-
bert W function. Suppose that z is real and positive. Then
√
Γ(z + 1) = (z + a(z))(z+1/2) e−(z+a(z)) 2π ,

so that
Γ(z + 1)
√ = (z + a(z))(z+1/2) e−(z+a(z)) .
2π
Raising both side to the power 1/(z+1/2) and then multiplying through
by −1/(z + 1/2) gives
1/(z+1/2)
−1 Γ(z + 1) z + a(z) − z+a(z)
√ =− e z+1/2
z + 1/2 2π z + 1/2
whence
1/(z+1/2) !
z + a(z) −1 Γ(z + 1)
− = Wk √
z + 1/2 z + 1/2 2π

and
1/(z+1/2) !
−1 Γ(z + 1)
a(z) = −z − (z + 1/2) Wk √ . (10.5)
z + 1/2 2π

137
Chapter 10. Consequences of Lanczos’ Paper

For real z the branches W−1 and W0 are real, yielding two solutions
which are here denoted a(−1, z) and a(0, z). Unfortunately equation (10.5)
is not of much practical use computationally without some additional
approximation of Γ itself, and the question of whether these equations
gives rise to more efficient Γ approximation schemes was not studied
further.
Equation (10.5) does, however, give a picture of what a(−1, z)
and a(0, z) look like, thanks to the Lambert W evaluation routines
of Maple 8. See Figure 10.1 for plots of these functions over the scaled
real line z/(1 − z), 0 ≤ z < 1. It is interesting to observe how little
these function appear to vary over the entire real line, in particular
a(−1, z) which is approximately 0.82 at z = 0 and decreases to about
0.79 as z → ∞.

a(−1,z/(1−z))
− 12 W−1 ( −1
π )

a(0,z/(1−z))
− 12 W0 ( −1
π )

0
0 z 1
z z
Figure 10.1: Path of a(0, 1−z ), a(−1, 1−z ), 0 ≤ z < 1

138
Chapter 11
Conclusion and Future Work

From the work done here it is clear that Lanczos’ formula merits se-
rious consideration as a viable alternative to Stirling’s series for the
efficient computation of the classical gamma function. Furthermore,
the generalized idea of the Lanczos transform described in Chapter 10
holds some promise as an approximation method for a larger class of
functions. There are, however, several unresolved issues and areas in
need of further investigation. This concluding chapter comments on
a number of unsettled questions, and conjecture answers to some of
these. Specifically:
With respect to Lanczos’ paper itself:

1. What is the precise asymptotic bound on the coefficients ak (r) of

the main series Sr (z) as k → ∞?
2. Can the convergence of Sr (z) be extended from Re(z) > −r to
Re(z) > −r − 1/2 as Lanczos claims?
3. How can the proof of the Lanczos Limit Formula (7.1) be extended
to Re(z) < 0?

With respect to the error bounding methods of Chapter 8:

1. Is there an analytical method to bound the relative error function

|ǫr,n (it/(1 − t))| for 0 ≤ t < 1, and hence |ǫr,n (z)| uniformly in
the right-half plane? More simply, in the n = 0 case, is there an
analytic bound for |ǫr,0 (it)| = |Fr (it) − 1|?
2. What is the precise relationship between the largest zero of the
relative error at infinity, ǫ∞
r,n , and the maximum of |ǫr,n (it/(1−t))|
on 0 ≤ t < 1?

139
Chapter 11. Conclusion and Future Work

3. How many zeros does ǫ∞r,n have? How are they distributed? What
is the largest one? The same questions apply to the individual
coefficients ak (r).

Finally, and perhaps most importantly, with respect to a deterministic

algorithm:

1. Is there a deterministic algorithm which, given z and ǫ > 0, pre-

scribes n and r so that Γ(z + 1) is computed with relative error
at most ǫ?

2. What is the relationship between r, n and the location of the

Lanczos shelf? How can the Lanczos shelf be used to best select
r as a function of n to minimize the relative error?

11.1 Outstanding Questions Related to Lanc-

zos’ Paper
Lanczos does not comment on the precise rate of decay of the coeffi-
cients ak (r) of the main series Sr (z). He does, however, make state-
ments about the rate and region of convergence of the series, both con-
sequences of the decay rate. The analysis of Chapter 4 determined
that the coefficients ak (r) = O (k −2r ) as k → ∞. Lanczos’ state-
ments suggest, and numerical evidence supports, the slightly faster rate
ak (r) = O (k −2r−1 ).
Closely related to this problem is whether convergence of the series
can be extended from the region Re(z) > −r found here to Re(z) >
−r−1/2 as claimed by Lanczos. If in fact ak (r) = O (k −2r−1 ) as k → ∞,
the proof of Theorem 5.1 would yield Lanczos’ bound. Lanczos’ claim
of convergence on the larger region seems so matter of fact that it
appears based on a simple principle, yet it is not clear what that is.
The analytic properties of the function represented by Sr (z), namely

Fr (z) = Γ(z + 1)(z + r + 1/2)−z−1/2 ez+r+1/2 (2π)−1/2 ,

suggests that in a manner similar to the Landau-Pringsheim Theo-

rem [1, p.237] for Dirichlet series, the abscissa of convergence can
be extended left until the first singularity (which is not a pole) at
z = −r − 1/2 is encountered. Perhaps a similar principle is at play

140
Chapter 11. Conclusion and Future Work

here. In the end, the asymptotic behaviour of the ak (r) and consequent
region of convergence of Sr (z) is a minor problem since, for computa-
tional purposes, only convergence on Re(z) ≥ 0 is required, and this is
guaranteed provided r ≥ 0.
Also of more theoretical interest than practical is Lanczos’ state-
ment of the limit formula of Chapter 7. The proof there established
the validity of the formula in the right-half plane Re(z) ≥ 0, while ac-
cording to Lanczos the formula is valid for all z ∈ C away from the
negative integers. His claim appears to be based on the asymptotic
behaviour of the coefficients ak (r) as r → ∞. Indeed, convergence of
the series Sr (z) extends further and further left in the complex plane
as r increases, but just how
∞
√ z+1/2 −(z+r+1/2)
X ′
lim 2π (z + r + 1/2) e ak (r)Hk (z)
r→∞
k=0
∞
X ′ 2 /r
= 2 lim r z (−1)k e−k Hk (z)
r→∞
k=0

directly is not entirely clear.

11.2 Outstanding Questions Related to Er-

ror Bounds
The main shortcoming of the Lanczos gamma approximation is the lack
of simple analytical error bounds on the resulting estimates. Instead,
for each non-negative integer n and r ≥ 0, the properties of ǫr,n (z) as an
analytic function of z on Re(z) ≥ 0 furnish simple empirical (uniform)
bounds which hold up when tested numerically. Nonetheless, to develop
an effective algorithm, one would like an a priori bound on ǫr,n (z) as a
function of its three parameters.
From the numerical evidence presented in Chapter 8 it is clear that
given n one should select r carefully in order to optimize the uniform
error bound. In all cases examined, setting r equal the largest zero of
ǫ∞
r,n results in uniform error bounds much improved on Lanczos’ own
estimates. Why is this so? It is not clear why this choice of r should
have such a dramatic effect on the uniform bound.

141
Chapter 11. Conclusion and Future Work

The nature of ǫr,n (z) as a function of r with z fixed (finite or infinite)

is not well understood, except in the simple n = 0 case of Section 10.3.
There we saw it is possible in the case of real positive z to select r which
makes ǫr,n (z) = 0, and consequently the approximation to Γ(z + 1)
becomes exact. Empirical evidence suggests that this is true for values
of n greater than zero as well. Refer to Figure 11.1 for a plot of ǫr,2 (z),
−0.5 < r < 3, −0.5 < z < 3 which shows several zeros of ǫr,2 (z) for
each value of z. For Re(z) > −1/2 and n ≥ 0 a fixed integer, if r(z)

0.4

−0.4
0 0

z r
3 3

Figure 11.1: ǫr,n (z), −0.5 < r < 3, −0.5 < z < 3

is defined to be a zero of ǫr,n (z), equation (10.5) and Figure 11.1 lead
one to conjecture that r(z) is a multivalued complex analytic function
with number of branches equal to the number of zeros of ǫ∞ r,n .

The data from Tables C.1 and C.2 suggests that neither the number
nor the distribution of ǫ∞r,n zeros follows a simple law. Even more per-
plexing is the striking similarity between the functions ǫ∞r,n and an+1 (r)
as illustrated in Figures 8.4, 8.5 and 8.6. This same behaviour was

142
Chapter 11. Conclusion and Future Work

found in all cases examined.

11.3 Conjectures on a Deterministic Algo-

rithm
In this section a deterministic algorithm is proposed for calculating
the gamma function using the truncated Lanczos formula. Here, a
deterministic algorithm is one which, given z with Re(z) ≥ 0 and ǫ > 0
as input, prescribes optimal n and r and computes Γ(1+z) with relative
error at most ǫ. As noted, the Lanczos method lacks a simple functional
relationship between the relative error ǫr,n (z) and the input parameters
n, r and z. Consequently, a deterministic algorithm is out of reach.
Even if the z parameter is removed from the requirement and we content
ourselves with a uniform error bound, the problem is still difficult.
In the uniform case, however, the empirical evidence strongly sug-
gests a simple deterministic algorithm. This evidence is the data of
Tables C.1 and C.2 of Appendix C, and the Lanczos shelf phenomenon
illustrated in Figures 1.5, 1.6 and 1.7 of Chapter 1. The similarity be-
tween Spouge’s and Lanczos’ methods noted in Section 2.6 would lead
one to guess that the error estimates ought to be similar as well. In
Spouge’s method, calculating Γ(1 + z) with uniform relative error ǫ > 0
requires O (− log ǫ) terms of the sum, as given by the bound (2.19). In
the Lanczos case, a plot of pairs (n, − log Mr(n),n ) from the tables of
Appendix C shows a near perfect linear relationship; see Figure 11.2.
A least squares fit of this data yields

n ≈ −2 − 0.3 log Mr(n),n .

Thus given ǫ > 0, one should choose n = ⌈−2 − 0.3 log ǫ⌉. The corre-
sponding r should then be chosen to be the largest zero of ǫ∞ r,n , which
is about size n. This choice of r ≈ n is consistent with shelf behaviour
demonstrated in Figures 1.5, 1.6 and 1.7. In each case, the transition
in decay rate of the ak (r) occurs at the shelf where r is approximately
equal to n.

143
Chapter 11. Conclusion and Future Work

210

− log Mr(n),n,15

6.8+3.3n

0
0 n 60

Figure 11.2: Linear fit of (n, − log Mr(n),n ), n = 0, . . . , 60

144
Bibliography
[1] T. Apostol. Introduction to Analytic Number Theory. Springer,
New York, 1998.

[2] E. Artin. The Gamma Function. Holt, Rinehard and Winston,

New York, 1964.

[3] J.P. Boyd. Chebyshev and Fourier Spectral Methods. Dover Pub-
lications, Inc., New York, 2000.

[4] F. Cajori. A History of Mathematical Notations. The Open Court

Publishing Company, Chicago, 1929.

[5] R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and

D. E. Knuth. On the Lambert W function. Adv. Comput. Math.,
5(4):329–359, 1996.

[6] P. J. Davis. Leonhard Euler’s integral: A historical profile of the

gamma function. The American Mathematical Monthly, 66:849–
869, 1959.

[7] Serge Dubuc. An approximation of the gamma function. J. Math.

Anal. Appl., 146(2):461–468, 1990.

[8] H.M. Edwards. Riemann’s Zeta Function. Academic Press, 1974.

[9] Paul Godfrey. A note on the computation of the conver-

gent lanczos complex gamma approximation. Web published at
http://my.fit.edu/~gabdo/paulbio.html, 2001.

[10] D. Gronau. Why is the gamma function so as it is? TEACHING

MATHEMATICS AND COMPUTER SCIENCE, 1:43–53, 2003.

[11] C. Lanczos. Trigonometric interpolation of empirical and analyt-

ical functions. Journal of Mathematics and Physics, 17:123–199,
1938.

[12] C. Lanczos. Applied Analysis. Prentice Hall, New Jersey, 1956.

[13] C. Lanczos. Linear Differential Operators. D. Van Nostrand Com-

pany, Ltd., New York, 1961.

145
Bibliography

[14] C. Lanczos. A precision approximation of the gamma function. J.

Soc. Indust. Appl. Math. Ser. B Numer. Anal., 1:86–96, 1964.

[15] Cornelius Lanczos. Collected published papers with commentaries.

Vol. I–VI. North Carolina State University, College of Physical
and Mathematical Sciences, Raleigh, NC, 1998. With a preface
by William R. Davis, a foreword by George Marx and a biograph-
ical essay by Barbara Gellai, Translated by József Illy, Laurent
Choquet, Don Ridgeway and Judith Kontság Meskó.

[16] S. Lang. Complex Analysis, Third Edition. Spring-Verlag, New

York, 1993.

[17] Y. Luke. The Special Functions and Their Approximations, Vol-

ume 1. Academic Press, New York, 1969.

[18] J.C. Mason and D.C. Handscomb. Chebyshev Polynomials. Chap-

man & Hall/CRC, Boca Raton, USA, 2003.

[19] P. Mikusinski and D. Taylor. An Introduction to Multivariable

Analysis from Vector to Manifold. Birkhäuser, Boston, USA, 2002.

[20] Edward W. Ng. A comparison of computational methods and

algorithms for the complex gamma function. ACM Trans. Math.
Software, 1(1):56–70, 1975.

[21] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling.
Numerical Recipes in C, Second Edition. Cambridge University
Press, Cambridge, UK, 1992.

[22] T.J. Rivlin. Chebyshev Polynomials from Approximation Theory

to Algebra and Number Theory. John Wiley & Sons, Inc., New
York, 1990.

[23] W. Rudin. Principles of Mathematical Analysis. McGraw-Hill,

Inc., New York, 1976.

[24] W. Rudin. Real and Complex Analysis. McGraw-Hill, Inc., New

York, 1987.

[25] L. Schwartz. Mathematics for the Physical Sciences.

Hermann/Addison-Wesley, Reading, Mass, 1966.

146
Bibliography

[26] Robert Spira. Calculation of the gamma function by Stirling’s

formula. Math. Comp., 25:317–322, 1971.

[27] John L. Spouge. Computation of the gamma, digamma, and

trigamma functions. SIAM J. Numer. Anal., 31(3):931–944, 1994.

[28] Nico M. Temme. Special Functions, An Introduction to the Classi-

cal Functions of Mathematical Physics. John Wiley & Sons, Inc.,
New York, 1996.

[29] S. C. van Veen. Review of a precision approximation of the gamma

function. Mathematical Reviews 31 #390, 1966.

[30] E.T. Whittaker and G.N. Watson. A Course of Modern Analy-

sis, Fourth Edition. Cambridge University Press, Cambridge, UK,
1965.

147
Appendix A
Preliminary Steps of Lanczos
Derivation

Reproduced here is the circuitous sequence of steps Lanczos’ uses to

arrive at (3.4) beginning with (3.1.1).
From Z ∞
z+1
Γ(z + 1) = α tz e−αt dt , Re(α) > 0 ,
0
make the substitution
α = 1 + ρz ,
where ρ > 0. This gives
Z ∞
z+1
Γ(z + 1) = (1 + ρz) (te−ρt )z e−t dt .
0

Next introduce the factor (eρ)−z

Z ∞
z+1 −z
Γ(z + 1) = (1 + ρz) (eρ) (ρte1−ρt )z et dt ,
0

and make the replacement v = e1−ρt to get

Z e
z+1 −z−1/ρ
Γ(z + 1) = (1/ρ + z) e [v(1 − log v)]z v 1/ρ−1 dv .
0

Finally, let r = 1/ρ − 1 to arrive at

Z e
z+1 −(z+r+1)
(z + r + 1) e [v(1 − log v)]z v r dv ,
0

which is equation (3.4).

148
Appendix B
A Primer on Chebyshev
Polynomials

The Chebyshev polynomials are a class of orthogonal polynomials (with

respect to a certain inner product) which are very effective for approx-
imation purposes. The literature on Chebyshev polynomials is vast, as
is the multitude of identities and relationships concerning these objects.
Stated here are the properties required in this work.
First, it should be noted there are several kinds of Chebyshev poly-
nomials, appropriately named Chebyshev polynomial of the first kind,
of the second kind, etc. The first two kinds are denoted {Tn (x)} and
{Un (x)}; we will be concerned only with this first kind, which is the
only polynomial solution for integer n of the differential equation

d2 y dy
(1 − x2 ) 2
− x + n2 y = 0 .
dx dx

The defining relationship for the {Tn (x)} is

Tn (x) = cos (nθ) where x = cos θ ,

and where x ranges from −1 to 1 as θ ranges from π to 0. The resulting

graph of Tn (x) oscillates between −1 to 1 and looks like a distorted
version of cos (nθ) . The leading coefficient of Tn (x) is 2n−1 and the
function is even or odd according as n is even or odd. The first few
Chebyshev polynomials of the first kind are:

T0 (x) = 1

T1 (x) = x

149
Appendix B. A Primer on Chebyshev Polynomials

T2 (x) = 2 x2 − 1
T3 (x) = 4 x3 − 3 x
T4 (x) = 8 x4 − 8 x2 + 1
T5 (x) = 16 x5 − 20 x3 + 5 x .

The family {Tn (x)}∞

n=0 is orthogonal with respect to the inner prod-
uct Z 1
dx
hf, gi = f (x)g(x) √ (B.1)
−1 1 − x2
and forms a basis for L2 [−1, 1] with norm induced by the inner product.
Further, if f (x) is differentiable on [−1, 1], (one sided differentiable at
the end points), then it is of bounded variation and may be written
∞
c0 (r) X
f (x) = + cn (r)Tn (x)
2 n=1

where convergence of the series is uniform. The coefficients cn (r) are

given by
2 1 dx
Z
cn (r) = f (x)Tn (x) √ . (B.2)
π −1 1 − x2

As an operational detail, take note that the family {Tn (x)}∞

n=0 is not
an orthonormal set, but only orthogonal. The proper normalization is
( r )∞
1 2
√ T0 , Tn (x) ,
π π
n=1

which must be taken into account for manipulations in the Hilbert space
setting.
The link between Chebyshev and Fourier series is the transformation
x → cos θ. This transforms f (x) into an even function of θ which can
then be written as a Fourier cosine series.
As noted, there are many identities involving Chebyshev polynomi-
als which follow from their interpretation as trigonometric functions.
The important ones for our purposes are:
√
T2n ( 1 − x2 ) = (−1)n T2n (x) (B.3)

150
Appendix B. A Primer on Chebyshev Polynomials

and the trigonometric form

T2n (sin θ) = (−1)n cos (2nθ) . (B.4)

Finally, the following important property of Tn (x) should be noted:

Of all monic polynomials of degree n on [−1, 1], the one with the least
absolute deviation from zero is Tn (x)/2n−1 . In this context, Tn (x)/2n−1
is the best approximation to the zero function on [−1, 1] by a monic
polynomial of degree n.
Refer to [17], [18] and [22] for proofs and a thorough treatment of
Chebyshev polynomials.

151
Appendix C
Empirical Error Data

Tables C.1 and C.2 summarize the results of numerical investigations

carried out on the relative error functions ǫr,n (z) and ǫ∞ r,n . For each
∞
n = 0, . . . , 60, the number of zeros of ǫr,n were counted, and the largest
and smallest zeros of were found to six decimals using a simple variation
of the bisection method. These statistics are listed in the first four
columns of the tables.
The last three columns of each table summarize data on the esti-
mated maximum of ǫr,n (it), 0 ≤ t < ∞ based on the theory of Sec-
tion 8.5. For r equal the largest zero of ǫ∞
r,n , the maximum of
∞
X
|ǫr,n (it)| = ak (r)Hk (it)


k=n+1

was estimated using the first 5 and then 15 terms of the sum. These
columns are labeled Mr(n),n,5 and Mr(n),n,15 , respectively. In addition,
Pn+15
the location of the maximum of
n+1 ak (r)Hk (it) was determined;

this data is in the last column labeled tmax .

152
Appendix C. Empirical Error Data

No. Smallest Largest

n Zeros Zero Zero r(n) Mr(n),n,5 Mr(n),n,15 tmax
0 2 −.223086 .319264 5.9 × 10−3 5.5 × 10−3 .88
1 4 −.173495 1.489194 1.0 × 10−4 1.0 × 10−4 2.13
2 6 −.151082 2.603209 6.4 × 10−7 6.3 × 10−7 4.27
3 8 −.137917 3.655180 8.5 × 10−8 8.5 × 10−8 4.05
4 8 −.129067 4.340882 4.3 × 10−9 4.3 × 10−9 5.05
5 10 −.122605 5.581000 1.2 × 10−10 1.2 × 10−10 6.34
6 12 −.117620 6.779506 2.7 × 10−12 2.7 × 10−12 7.88
7 14 −.113619 7.879012 3.9 × 10−14 3.9 × 10−14 6.72
8 14 −.110313 8.406094 6.9 × 10−15 6.9 × 10−15 9.23
9 16 −.107519 9.656578 2.1 × 10−16 2.1 × 10−16 10.51
10 18 −.105114 10.900511 6.1 × 10−18 6.1 × 10−18 11.83
11 20 −.103013 12.066012 1.1 × 10−19 1.1 × 10−19 14.30
12 22 −.101157 13.144565 5.2 × 10−21 5.2 × 10−21 12.38
13 22 −.099499 13.726821 4.0 × 10−22 4.0 × 10−22 14.69
14 24 −.098005 14.977863 1.2 × 10−23 1.2 × 10−23 15.98
15 26 −.096650 16.209805 3.6 × 10−25 3.6 × 10−25 17.45
16 28 −.095412 17.345444 3.1 × 10−27 3.1 × 10−27 23.47
17 30 −.094275 18.399283 5.0 × 10−28 5.0 × 10−28 18.11
18 30 −.093226 19.048512 2.5 × 10−29 2.5 × 10−29 20.16
19 32 −.092252 20.298892 7.8 × 10−31 7.8 × 10−31 21.45
20 34 −.091346 21.508926 2.1 × 10−32 2.1 × 10−32 23.40
21 36 −.090499 22.618910 1.8 × 10−34 1.8 × 10−34 17.14
22 36 −.089704 23.118012 5.2 × 10−35 5.2 × 10−35 24.34
23 38 −.088958 24.370498 1.7 × 10−36 1.7 × 10−36 25.62
24 40 −.088253 25.617904 5.2 × 10−38 5.2 × 10−38 26.97
25 42 −.087588 26.798597 1.1 × 10−39 1.1 × 10−39 30.13
26 44 −.086957 27.886311 3.6 × 10−41 3.6 × 10−41 24.97
27 44 −.086358 28.440357 3.5 × 10−42 3.5 × 10−42 29.80
28 46 −.085789 29.692534 1.1 × 10−43 1.1 × 10−43 31.09
29 48 −.085246 30.931341 3.4 × 10−45 3.4 × 10−45 32.63
30 50 −.084727 32.080670 4.4 × 10−47 4.4 × 10−47 39.41

Table C.1: Error Data, n = 0, . . . , 30

153
Appendix C. Empirical Error Data

No. Smallest Largest

n Zeros Zero Zero r(n) Mr(n),n,5 Mr(n),n,15 tmax
31 52 −.084232 33.145772 4.2 × 10−48 4.2 × 10−48 31.63
32 52 −.083757 33.762726 2.4 × 10−49 2.4 × 10−49 35.27
33 54 −.083302 35.014250 7.7 × 10−51 7.7 × 10−51 36.57
34 56 −.082865 36.235367 2.2 × 10−52 2.2 × 10−52 38.65
35 58 −.082445 37.356480 7.5 × 10−55 7.5 × 10−55 34.71
36 60 −.082041 38.385241 3.8 × 10−55 3.8 × 10−55 38.30
37 60 −.081651 39.085095 1.7 × 10−56 1.7 × 10−56 40.74
38 62 −.081275 40.334630 5.3 × 10−58 5.3 × 10−58 42.09
39 64 −.080912 41.529155 1.3 × 10−59 1.3 × 10−59 45.44
40 66 −.080562 42.626437 2.7 × 10−61 2.7 × 10−61 36.11
41 66 −.080223 43.154830 3.6 × 10−62 3.6 × 10−62 44.92
42 68 −.079894 44.407411 1.2 × 10−63 1.2 × 10−63 46.21
43 70 −.079576 45.651117 3.7 × 10−65 3.7 × 10−65 47.74
44 72 −.079268 46.814382 6.2 × 10−67 6.2 × 10−67 54.32
45 74 −.078969 47.889652 3.8 × 10−68 3.8 × 10−68 44.46
46 74 −.078679 48.477371 2.6 × 10−69 2.6 × 10−69 50.39
47 76 −.078396 49.729491 8.2 × 10−71 8.2 × 10−71 51.68
48 78 −.078122 50.959691 2.5 × 10−72 2.5 × 10−72 53.73
49 80 −.077855 52.092791 1.7 × 10−74 1.7 × 10−74 71.81
50 82 −.077596 53.141340 3.8 × 10−75 3.8 × 10−75 51.74
51 82 −.077343 53.799879 1.8 × 10−76 1.8 × 10−76 55.85
52 84 −.077096 55.050733 5.8 × 10−78 5.8 × 10−78 57.19
53 86 −.076856 56.257932 1.5 × 10−79 1.5 × 10−79 60.47
54 88 −.076622 57.365268 1.8 × 10−81 1.8 × 10−81 44.65
55 88 −.076393 57.869538 4.0 × 10−82 4.0 × 10−82 60.03
56 90 −.076170 59.122331 1.3 × 10−83 1.3 × 10−83 61.32
57 92 −.075952 60.369399 4.1 × 10−85 4.1 × 10−85 62.81
58 94 −.075739 61.546699 8.5 × 10−87 8.5 × 10−87 68.92
59 96 −.075531 62.631604 3.3 × 10−88 3.3 × 10−88 56.29
60 96 −.075327 63.192152 2.9 × 10−89 2.9 × 10−89 65.50

Table C.2: Error Data, n = 31, . . . , 60

154