UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

An analysis of the Lanczos gamma approximation Pugh, Glendon Ralph 2004

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2005-995364.pdf [ 8.16MB ]
Metadata
JSON: 831-1.0080001.json
JSON-LD: 831-1.0080001-ld.json
RDF/XML (Pretty): 831-1.0080001-rdf.xml
RDF/JSON: 831-1.0080001-rdf.json
Turtle: 831-1.0080001-turtle.txt
N-Triples: 831-1.0080001-rdf-ntriples.txt
Original Record: 831-1.0080001-source.json
Full Text
831-1.0080001-fulltext.txt
Citation
831-1.0080001.ris

Full Text

A N A N A L Y S I S O F T H E L A N C Z O S G A M M A A P P R O X I M A T I O N by G L E N D O N R A L P H P U G H B.Sc. (Mathematics) University of New Brunswick, 1992 M.Sc . (Mathematics) University of Br i t i sh Columbia, 1999 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S (Department of Mathematics) We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A November 2004 © Glendon Ralph Pugh, 2004 Abstract This thesis is an analysis of C . Lanczos' approximation of the classical gamma function T(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 proof's of al l 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. A t 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 s t i l l widely used today. Another more recent method is that of Spouge [27], which is similar in form though different in origin than Lanczos' formula. A l l three of these approximation methods take the form ; r(z + l) = V2^(z + wy+V2e-z-w[sw<n(z) + eWtn(z)} (1) where sw>n(z) denotes a series of n + 1 terms and ew,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, wi th w > 0 fixed, as n —> oo the series sw,n(z) converges while eWtn(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 | e W ) n ( z ) | in the right half plane based on the behaviour of this function as \z\ —> oo. 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 I _ p - l 2 / r _ L _ + r - 2 3 / r * ( * ~ ! ) , T(z + 1) = 2 l im r r—>oo as stated without proof by Lanczos at the end of his paper. i i Table of Contents Abst rac t i i Table of Contents i i i L is t of Tables v i L is t of Figures v i i Acknowledgments v i i i Chapter 1. In t roduct ion 1 1.1 Lanczos and His Formula 2 1.2 Thesis Overview 8 Chapter 2. A P r i m e r on the G a m m a Funct ion 20 2.1 First Encounters wi th Gamma 20 2.2 A Brief History of Gamma 22 2.3 Definition of T 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 Addi t iona l 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{6) 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 Addi t iona l 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 i i i Table of Contents 4.1.2 Closed Form for fr(x) and fE,r{x) 54 4.2 Smoothness of vr(x), fr(x) and fE,r{x) 56 4.3 Smoothness of vr(9), fr(0) and fE]r(9) 59 4.3.1 Generalized fr(0) 59 4.3.2 Derivatives of G(0, j,k,l;d) 61 4.3.3 Appl icat ion to fr 61 4.4 Fourier Series and Growth Order of Coefficients 64 Chapter 5. A n a l y t i c Cont inuat ion of M a i n Formula 66 5.1 Domain of Convergence of Sr(z) 66 5.2 Analy t ic Continuation 69 5.3 Addi t ional Remarks 69 Chapte r 6. C o m p u t i n g 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 Par t ia l Fraction Decomposition 80 6.7 Ma t r i x Representation 82 Chapter 7. Lanczos ' L i m i t Formula 84 7.1 Some Motivat ion 85 7.2 Proof of the Limi t Formula 87 7.3 Addi t ional Remarks 95 Chapte r 8. Implementat ion & E r r o r 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 e™n 109 8.5.1 The Motivat ing Idea .' 109 8.5.2 Zeros of e^n I l l iv Table of Contents 8.5.3 Improved Lanczos Formula ' I l l 8.5.4 Opt imal Formulas 114 8.6 Addi t ional Remarks 116 Chapter 9. Compar ison of the Methods 120 9.1 Stirling's Series 120 9.1.1 r(20 + 17i) wi th \ep\ < l O " 3 2 120 9.1.2 T(l + z) wi th \ep\ < 10~ 3 2 Uniformly 122 9.2 Spouge's Method 123 9.2.1 T(20 + 17i) wi th \ep\ < 10~ 3 2 123 9.2.2 T ( l + z) wi th \ep\ < 10~ 3 2 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. Conclus ion and Future W o r k 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 Algor i thm 143 Bib l iography 145 A p p e n d i x A . P re l imina ry Steps of Lanczos Der iva t ion 148 A p p e n d i x B . A P r i m e r on Chebyshev Polynomials 149 A p p e n d i x C . E m p i r i c a l E r r o r D a t a 152 List of Tables 1.1 Relative decrease of dk(r), r = 1,4, 7 4 2.1 Upper bound on |i?jv,n(6 + 13i)| in Stirl ing Series 33 2.2 Coefficients of Spouge's Series, a = 12.5 38 6.1 Scaled summands of c 6(6) 75 6.2 Coefficients as functions of r 80 8.1 Lanczos' uniform error bounds 106 8.2 Location of maximum of r]r,n{6) 107 8.3 Errors at infinity 107 8.4 Approximate Zeros of e~6 113 8.5 Sample coefficients, n = 1 0 , r = 10.900511 116 8.6 Comparison of M r ( n ) n and an+i(r(n)) 119 9.1 Upper bound on |£jv,n(19 + 17z)| in Stirl ing Series 121 9.2 Upper bound on |.Ejv,ra(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 . l Error Data, n = 0 , . . . , 30 153 C.2 Error Data, n = 3 1 , . . . , 60 154 v i List of Figures 1.1 Cornelius Lanczos •. 3 1.2 Relative decrease of a/»(l) 5 1.3 Relative decrease of ajk(4) : 5 1.4 Relative decrease of a,k(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{0), - T T / 2 <0<IX/2 13 2.1 r(> + 1), - 4 < z < 4 27 2.2 \T{x + iy)\, - 3 < x < 4, - 3 < y < 3 28 3.1 Pa th of v(l - logv) 44 4.1 v(x)/e, - 1 < x < 1 54 4.2 v / 2 / £ ; , r ( ^ ) / e r , - 1 < x < 1, r = 0 ,1 ,2 ,4 55 5.1 Growth of ak(r), k'^^, and k~W wi th r = 1.25 71 7.1 exp [-r{0 - TT/2) 2] - gr(0) for r = 1, 6,10 86 7.2 logQ(fc ,10), fc = 0 , . . . , 2 4 ' 9 6 7.3 logQ(fc ,20) , k = 0 , . . . , 2 4 97 7.4 logQ(/c,50), fc = 0 , . . . , 2 4 97 8.1 [ e ~ ] 1 / 5 , n = 0 , l , 2 , 3 , - l / 2 < r < 6 112 8.2 log | e P , ) 6 ( i t / ( l - * ) ) ]> 0 < t < l , j = 0 , . . . , 1 1 114 8.3 log \e r A ( i t /{ l - t))\, 0 < t < 1, 0 < r < 5 115 8.4 [e^]1/'?, - 0 . 2 < r < 7 117 8.5 M 0 ] 1 / 7 , -0 -2 < r < 7 118 8.6 [e~6]x/7 and M r ) ] 1 / 7 , - 0 . 2 < r < 7 118 10.1 Pa th of a(0,^), a ( - l , I f I ) , 0 < z < 1 138 11.1 e r , n (z), - 0 . 5 < r < 3, - 0 . 5 < z < 3 , 1 4 2 11.2 Linear fit of (n, - log M r ( n ) , n ) , n = 0 , . . . , 60 144 vi i Acknowledgments I would like to thank al l 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 . B i l l Casselman for his patience, guidance and encouragement, my committee members Dr . David Boyd and Dr . Richard Froese, and Ms . Lee Yupi tun 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 Nor th Carol ina State Un i -versity for his k ind 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, M P J A V A , a multiple precision floating point computation package in Java produced by The H A R P O O N Project at the University of Nor th Carol ina was used for many of the high precision numerical investigations. In addition, the online spreadsheet MathSheet developed by Dr . B i l l Casselman and Dr . David Aus t in 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. v i i i Chapter 1 Introduction The focus of this study is the little known Lanczos approximation [14] for computing the gamma function T(z + 1) wi th complex argument z. The aim is to address the following questions: (i) Wha t 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 T(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. W i t h i n a year of Euler 's first consideration of the gamma function in 1729, Stir l ing pub-lished the asymptotic formula for the factorial wi th 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 wi th optimal computa-tional strategies. There are, however, other methods for computing the gamma function wi th 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 Stirl ing 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 T(z-\-1) on Re(z) > 0. In the specific case of a uniformly bounded relative error of 1 0 - 3 2 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 wi th 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\ —» oo. (ii) A . careful examination of Lanczos' paper [14], including a proof of the l imit 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 T(z + 1) = v 7 ^ (z + r + l / 2 ) 2 + 1 / 2 e~(z+r+lW Sr(z) (1.1) where Sr(z) = 1 . . , . z , . z(z — 1) -ao(r) + fll(r) — + ^(0 ( z + 1 ) ( g + 2 ) + (1.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, i n a manner similar to a divergent series, the coefficients afc(r) ini t ial ly decrease quite rapidly, followed by slower decay as k increases. Furthermore, wi th 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 afc(r) i n (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-lln fact, 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 afc(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 \dk+i{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 | a f c + i ( l ) / a f e ( l ) | | a f c + 1 (4) /d f c (4) | |a f c + 1 (7) /a f c (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 a^(r), r = 1,4, 7 coefficients drops very quickly at k •= 1, but then flattens out into a more regular pattern. Compare this wi th 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+i(r)/a,k(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 ini t ial large decrease, the ratios do not follow a smooth pattern as in the r = 1 case, but rather, jump around irregu-larly. The afc(r) are actually Fourier coefficients of a certain function, and are bounded asymptotically by a constant times k~2r. However, the in i -4 Chapter 1. Introduction i i i 1 t ; ; 1 M M ! ! 1 i ! 1 1 1 1 1 1 1 I I ] ! i M ; ! ; j \ i i i i M i 1 I — r r r >6.la ~ j " " " ; T j -r j ' i i M i . j |! 1. , i . fcj-i/l)/^^)! , _ ; 1 ! ; I i "| ; 1 i l l : i t ; ' ! j • | | : f-Figure 1.2: Relative decrease of <2fc(l) 12, . ; |-l9g|aAii(4)/a f c(ji)J 29 Figure 1.3: Relative decrease of a^(4) t ia l decrease of the coefficients appears more exponential in nature, and then switches abruptly to behave according to the asymptotic bound. 5 Chapter 1. Introduction This peculiar phenomenon is nicely illustrated in Figures 1.5, 1.6 and 1.7 where — log \aic(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 /c-term approx-imation should use r ~ k, since the decay rate of coefficients after the cutoff point slows considerably. Apar t 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. A t 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 i— "-log|a f c(20)/e 2 0 | 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 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 al l state-8 Chapter 1. Introduction merits. (ii) From a practical point of view, to determine how (1.1) is best used to compute T(z +1), which reduces to an examination of the error resulting from truncation of the series (1.2), and to compare Lanczos' method wi th 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 T(z + 1): Chapter 2 II. Lanczos' Paper: Details and Proofs; Lanczos L imi t Formula: Chapters 3-7 III. Error Discussion: Chapter 8 I V . Comparison of Calculation Methods: Chapter 9 V . Consequences and Extensions of the Theory: Chapter 10 V I . 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, wi th the aim of providing the reader wi th a survey of the main ideas, thus serving as a guide for navigating the rest of the thesis. H i s t o r y and Background of T(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 9 Chapter 1. Introduction which is the form wi th which Lanczos begins his study. Following the formal definition of T(z'+ 1), a number of standard identities are given wi th some explanation but without formal proof. The chapter concludes wi th a discussion of computational methods, beginning wi th Stirling's series, which is a generalization of the familiar Stirling's formula n\ ~ \j2im nne n as n —> oo . Stirling's asymptotic series forms the basis for most computational al-gorithms for the gamma function and much has been written on its implementation. A n example on the use of Stirling's series is given wi th 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 T(z + !) = (* + a)z+1/2e^z+a\2TT)1/2 t 1 z + k k=l where N = \a] — 1, c 0 = 1 and ck is the residue of T(z + l)(z + a)-(z+i/2) ez+a^2 7 r)- 1/2 a t z — Al though 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 Cfc are much easier to compute. This method also extends to the digamma function ty(z + 1) = d/dz [logr(z + 1)] and trigamma function ^f'(z). Lanczos' Paper: Details and Proofs Chapter 3 examines in detail Lanczos' paper, beginning wi th the deriva-t ion 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 wi th a free parameter r > 0, and then replacing that function wi th its Fourier series. Specifically, beginning wi th poo F{z + 1) = / t'e-'dt , . , Jo 10 Chapter 1. Introduction these steps yield T(z + 1/2) = (z + r + l/2)z+1/2e^z+r+1^V2 rn/2 x / c o s 2 z 0 ' - T T / 2 V2v(6)rsm6 ; i .3) de \ogv(0) The function v(9) in the integrand is defined implici t ly by v(l — log v) = cos 2 9 where v is increasing, and v = 0,1, e corresponds to 6> = —TT/2, 0,7r/2, respectively. Denoting by fr(0) the term in square brackets in the integrand of (1.3), and JE,r(9) its even part, the properties of this last function are such that it can be replaced by a uniformly convergent Fourier series 1 fe=l Integrating this series against cos 2 2 9 term by term then gives rise to the series Sr(z) of (1.2) wi th 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 afc(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 JE,T{9) governs the rate of asymptotic decrease of the coefficients a^(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{9), fr(0) and JE,r{9) 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 a^(r) is established, a bound on the functions i 1 if k = 0, , x (z+l)---(z+fc) — 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 afc(r). The standard method provided by Fourier theory, namely the direct integration 2 r / 2 ak(r) = - fEr(6) cos (2k0) dO (1.5) proves impractical due to the complicated nature of fE,r(&)- Lanczos overcomes this difficulty wi th the clever use of Chebyshev polynomials as follows: define Tr{z) = 2-1/2r(z + l/2)(z + r + l/2)-z-1/2exp(z + r + l/2) , which is just the integral in (1.3), and recall the defining property of the 2kth Chebyshev polynomial: k r2fe(cOS 6) = C2j,2k cos2j 0 3=0 = cos(2k0) . Then (1.5) becomes ak(r) = - f ^ fE,M E C ^ co s2 j 9 d 0 2 ^ = -^C2j,2kfr{j) , 7T * — ' 3=0 a simple linear combination of Tr 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 wi th the following formula which he claims is the l imit ing behaviour of (1.1), and which defines the gamma function in the entire complex plane away from the negative integers: T{z + 1) = 2 l im rz r—>oo -. oo ^ + E(-1)/£e_fc2/r^) fc=l ; i .6) 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 l imit (1.6) is proved in detail, a result termed the "Lanczos L imi t Formula" in this work. The proof is motivated by the plots in Figure 1.8 of the function fr(0) from equation (1.3). Notice that as r increases, the area un-, % / 2 e 2 r r = 2 r = l r=0 Figure 1.8: fr(0), -TT/2 <6<ir/2 13 Chapter 1. Introduction der fr(0) appears to be increasingly concentrated near 9 = n/2 where / P ( T T / 2 ) = V2er. _ The l imit formula (1.6) suggests a rescaling of fr(6) to y/r/(2iT)fer(6)e which decays rapidly to zero on [—7r/2,7r/2] except at 0 = ir/2 where it has the value y/r/n. This behaviour is reminiscent of that of the Gaussian distribution y/r/nexp [—r(6 — 7r/2)2] on (—oo,7r/2], and in-deed, plots of these two functions in the vicinity of 9 = ir/2 show a remarkable fit, even for small values of r . W i t h some work, the l imit ing value of (1.3) is shown to be expressible in the form /oo cos 2 2 d^e-r^9-n/2)2 dO -co which, upon replacement of cos 2 2 9 by its Fourier series, gives rise to the coefficients v ia cos (2k9)r1'2e-r{e-v,2)2 dt = ( - l ) V f c 2 / r . A comparison of (1.6) with the main formula (1.1) suggests that for large r, ,/f ~ < - l ) V - > , (1.7) and plots of these coefficients for various r values supports this rela-tionship. From the smoothness analysis of fr(9) it follows that a/-(r) = O (fc_r2rl). 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 wi l l be the value of the coefficients at which the convergence begins to slow down. A t the same time, however, we have to wait longer, before the asymptotic stage is reached. Al though 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 L imi t Formula does not clarify quantitatively Lanczos' observation, but it does shed some light on what perhaps motivated h im 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 wi th 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\ —> oo. 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 wi th Re(z) > 0 since the reflection formula provides estimates for gamma in the left hand plane if the function in right half plane is known. For notation, write Sr(z) = Srtn(z) + £r,n(z) where Sr^n(z) is the sum of the first n terms and Following Lanczos, the function er^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 |e 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 |e r > n (z) | , observe that ertn(z) is an analytic func-t ion on Re(z) > 0. It turns out that oo fc=n+l 15 Chapter 1. Introduction a constant denoted e^ , the error at infinity. Under these conditions, as a consequence of the maximum modulus principle, the maximum of l er,n(z)| o n Re(z) > 0 must occur on the line Re(z) = 0, possibly at infinity. Furthermore, thanks to the Schwartz reflection principle, this maximum wi l l occur on the positive imaginary axis. Thus to bound l er,n(z)| o n Re(z) > 0, it is sufficient to bound | e r ) n ( i£ ) | , 0 < t < oo: st i l l a daunting task.. The problem is overcome by first transforming the domain to a finite interval v ia the mapping z(t) = it/(l — t), 0 < t < 1, and then observing that under this mapping the functions Hk(z(t)) have the simple well behaved form Thus the maximum of |e r ) n(,z)| on Re(z) > 0 can be easily estimated empirically by examining the first few terms of \erin(z(t)) \ on 0 < t < 1. This estimation method is used to examine the dependence of er,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, wi th 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 e^°n = 0 have solutions? This question motivated an extensive numerical examination of the relative error functions er^n(z) and e£°n for n = 0 , . . . , 60. The results of this investigation are tabulated in Appendix C . Experimentally, e^n = 0 was found to have a finite number of real solutions. For the values 0 < n < 60 examined, e^°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 was found to give the least uniform bound on |e 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 l e r , n ( ^ ) | < 1 1 0.001 1 1.5 0.00024 2 2 5.1 x 10~ 5 3 2 1.5 x 10~ 6 3 3 1.4 x 10~ 6 4 4 5 x 10~ 8 6 5 2 x 1 0 " 1 0 B y contrast, selecting r to be the largest zero of e^°n yields the following dramatically improved uniform error bounds (again, for Re(z) > 0): n r \€r,n{Z)\ < 1 1.489194 1.0 x 10" 4 2 2.603209 6.3 x 10" 7 3 3.655180 8.5 x 10" 8 4 4.340882 4.3 x 1 0 - 9 6 6.779506 2.7 x 10~ 1 2 Compar i son of Calcu la t ion Methods The purpose of this chapter is to compare the methods of Lanczos, Spouge and Stirl ing in the context of an extended precision compu-tation. For each method, a detailed calculation of T(20 + I7i) wi th relative error |e p | < 10~ 3 2 is carried out. Additionally, formulas for T ( l + z) wi th uniform error bound of 10~ 3 2 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 wi th 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(2) > 0 as F{z) = / cos2* 0o-(0) dO where g{6) e L2[—7r/2, TT/2] is even, then r^Oz + 1/2) ri z z(z-l) F i-Z ) = ^T(7TTT k ° + "l7+i+a*(z + 1 ) ( » + 2) + • ' ' The dk in the series are the Fourier coefficients of g(9), and these are given by 2 ^ ak = -J2 C2j,2kF(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 L imi t 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 T{z + 1) = (z + r + i / 2 ) ( z + 1 / 2 ) e - ( 2 + r + i / 2 ) ^ 2 T J (1 + €rfi{z)) , where r = 0.319264, the largest zero of e£°0, and |e r )o(^)| < 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 e r(2) io(^) = 0, thus yielding r(z + 1) = (z + r(z}+ l / 2 ) ( z + 1 / 2 ) e - ( z + r ( 2 ) + 1 / 2 ) v / 2 ^ : ? The answer is shown to be yes, wi th 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 wi th 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. Lett ing r(n) denote the largest zero of e^°n and Mr(n^n the maximum of |er,n(^)| f ° r 0 < t < oo, a plot of (n, — l o g M r ( n ) i n ) reveals a near perfect linear trend n = — a l o g M r ( n ) n + b. The algorithm is then clear: given e > 0, select n = \—aloge + 6], and set r = r(n). Then (1.1) truncated after the n t h term computes T ( l + z) wi th a uniformly bounded relative error of at most e. Al though 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 wi th a survey of two computational methods wi th 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 wi th 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 wi th 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 ini t ial exposure is accidental: n! is that button on their scientific calculator which causes overflow errors for al l 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 — 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 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 /•oo n ! = / fe^dt (2.1) Jo 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! ~ V2wnnne~n as n —> oo . A t 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 wi th the gamma function begin wi th a statement of the definition, usually in terms of an integral or a l imit . 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 + l) 1 + 2 +••• + n = V 7 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 unti l 1808 when C . K r a m p 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 + l)x xl = l im j-1 , 2.2 •m-»oo \ X + 1) • • • [X + m) and the second an integral representation: x\= f ( - l o g t f dt . (2.3) Jo 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= — logt puts (2.3) into the more familiar form In 1808 Legendre introduced the notation which would give the gamma function its name 1 . He defined so that r ( 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 T(n + 1) in-stead of r (n) is due to Legendre and void of any rationality. This unfortunate circumstance compels us to utilize the no-tation z\ instead of T(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 U(s) = s!, stating: Unfortunately, Legendre subsequently introduced the no-tation r(s) for n(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 T(x) as the Mel l in transform of exp (—t), so that Legendre's normalization is the right one in some Whittaker and Watson [30] give 1814 as the date, however. x\ = du . (2.4) 23 Chapter 2. A Primer on the Gamma Function sense, though Mel l in was not born unti l some twenty one years after Legendre's death. A n interesting recent investigation into the history of the T 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 > —I (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 T(x) by the Bohr-Mollerup theorem (see [23]): T h e o r e m 2 .1 . Suppose f : (0, oo) —> (0, oo) is such that / ( l ) = 1, f(x + 1) = xf(x), and f is log-convex. Then f = F. So as a function of a real variable, T(x) has nice behaviour. Bu t if gamma is extended beyond the real line, its properties are even more satisfying: as a function of a complex variable z, F(z) = / 0°° i 2 _ 1 e _ t dt defines an analytic function on Ke(z) > 0. Along wi th the fundamental recursion F(z + 1) = zF(z), F can then be extended to an analytic function on z £ C \ {0, —1, —2, —3, . . .} . Indeed, Euler's formula (2.2) converges and defines F 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 wi th certainty that a particular definition is in some way the correct one. In the present case, £l definition of the correct generalization of n! is sought. The Bohr-Mollerup theorem and the analytic properties of Euler's F strongly suggest that (2.4) is indeed the right definition. The appearance of Euler's F 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 al l kinds of problems. Moreover, the function has a number of interesting proper-ties. So is F(z) = J 0°° tz~le~l dt the natural extension of nl ? It seems so, and we wi l l 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 T A s 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] wi l l be used: D e f i n i t i o n 2.1. For z G C wi th Re(z) > —1, the gamma function is defined by the integral From this definition, we immediately state T h e o r e m 2.2. For T(z + 1) defined by (2.5): (i) T(z + 1) is analytic on Re(z) > —1; (ii) T(z + 1) = zY{z) on R e ( » > 0 ; (iii) T(z + 1) extends to an analytic function on C with simple poles at the negative integers. P r o o f o f T h e o r e m 2.2: (i) Let Q denote the domain Re(z) > —1. The integrand of tze~l dt is continuous as a function of t and z, and for fixed t > 0 is ana-lytic on f2. Further the integral converges uniformly on compact subsets of f2. Hence T(z + 1) is analytic on £1. (ii) For Re(z) > 0, integrating by parts gives (iii) For Re(z) > 0 and any integer n > 0, the fundamental recursion T(z + 1) = zT(z) implies (2.5) T(z) T(z + n+l) nz=o(*+*o' (2.6) 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 Y(z) on Re(z) > 0. Hence (2.6) is the unique analytic continuation of T(z) to the right half plane Re(z). > — n—1, z not a negative integer. Since n was an arbitrary non-negative integer, V(z) can be extended to an analytic function on C \ { 0 , - 1 , - 2 , - 3 , . . . } . For z = —n a negative integer, (2.6) shows that Y has a simple pole wi th residue . . Y(z + n +1) 2 h m n ( z + n)r(z) = hmjz + n) Y^fJTV) r(i) m:J(-^  + ^ ) 1 re:01(-^+^) (-ir _ A n interesting generalization of (2.5) due to Cauchy leads to another analytic continuation of the gamma function. Let o = . Re(—z — 1). Then ro + i)= r**fe_t- E • <2-7) A n empty sum in the integrand of (2.7) is to be considered zero. This integral converges and is analytic for z G C such that Re(z) ^ - 1 , - 2 , - 3 , . . . and reduces to (2.5) if Re(z) > - 1 . See [30, pp.243-244] Figure 2.1 shows a graph of Y(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 \Y(x + iy)\ which illustrates the behaviour of the function away from the real axis. 26 Chapter 2. A Primer on the Gamma Function Figure 2.1: T(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 T(z + 1) = zY(z) holds for al l z E 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 T{z) = l im — — (2.8) V ' n-^ oo Z(Z + 1) • • • (Z + n) V ' which converges for all z E C \ { 0 , — 1, — 2, — 3 , . . . } . This identity follows from integrating by parts in ( l — | ) n ^ 2 _ 1 dt to get n J z(z + 1) • • • (z + n) 27 Chapter 2. A Primer on the Gamma Function Figure 2.2: \F(x + iy)\, - 3 < x < 4, - 3 < y < 3 and then letting n —> oo . Rewrit ing (2.8) in reciprocal form as l im ^ + 1 ) ' / ' ( " + n ) = l im * e * - k « * * n - i i/*> ff f 1 + ^ e - / f c fc=l yields the Weierstrass product The 7 appearing in (2.9) is the Euler (or sometimes Euler- Mascheroni) constant 7 = l i m ^ ^ YJl=\ (V f c) _ l o S n ~ 0.5772 Wri t ing both T(z) and T(-z) using (2.9) and recalling the infinite 28 Chapter 2. A Primer on the Gamma Function product representation of the sine function s i n z — z n ( 1 _ fc%2 j fc=i ^ y we arrive at the very useful reflection formula r(i + z)r(l - z) = 7VZ (2.10) smnz 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 T(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 log r(z + i)H(z + k) {z + N + l/2)hg(z + N) k=i (2 .H) 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 AT = 0. The B2j are the Bernoulli numbers, which are defined by writ ing t/(et — 1) as a Maclaur in series: t oo R 3=0 J ' The function B2n(x) is the 2nth Bernoulli polynomial defined over K. by periodic extension of its values on [0,1]. B y judiciously selecting n and N the absolute error in logr(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 wi th the relative error of T(z + 1) itself. Equation (2.11) is valid in the so-called slit plane C \ {t e R 11 < -N}. Equation (2.11) wi th N = 0 is derived using Euler-Maclaurin sum-mation. Alternatively, again wi th N = 0, begin wi th the formula of Binet obtained from an application of Plana's theorem, logrOz + 1) = (-2 + 1/2) logz 1 r z + - log 2?r + 2 / 2 Jo arctan (t/z) dt and expand arctan (t/z) in the last term into a Taylor polynomial wi th 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) wi th N > 1 is obtained by replacing z wi th z + N and . applying the fundamental recursion. Exponentiation of Stirling's series (with N = 0) yields Stirling's asymptotic formula for T(z + 1) itself: T(z+l) = e-zzz+1/2(2ir 11/2 1 1 1 12z + 288z 2 + (2.12) 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: r(* + i) V2^zze~z as \z\ CO (2.13) It may be tempting to let n —> oo in (2.11) wi th the hope that the integral error term wi l l go to zero, but unfortunately the series contain-ing the B2j becomes divergent. The terms of this sum init ial ly decrease 30 Chapter 2. A Primer on the Gamma Function i n modulus, but then begin to grow without bound as a result of the rapid growth of the Bernoulli numbers. The odd Bernoulli numbers are al l zero except for B\ = —1/2, while for the even ones, Bo = 1, and for 3 ~ ^ E - (-l)J+12(2j)!C(2j) ( 2 u ) BV - ^ • (2.14) Since Q{2j) = }ZT=i k~2j -» 1 rapidly with increasing j , B2j = O ((2J)!(2TT) (see [8, p.105]). Wha t is true, however, is that the integral error term in (2.11) tends to zero uniformly as \z\ —> oo in any sector | arg (z + N)\ < S < TT. In terms of T(z + 1) itself, this means the relative error tends to zero uniformly as \z\ —• oo in any sector | arg (z + N)\ < S < 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]: T h e o r e m 2.3. Let 9 = axg(z + N), where —TT < 9 < TT, and denote the absolute error by J o 2n(z + N + x)2n Then B2n+2 . \ 2n+2 \EN,n(z)\ < ' ^ COS (0/2)) (2n + 2) (2n+ l)(z + N)2n+1 (2.15) 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 (9/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 E x a m p l e 2 .1 . Compute r(7+13i) accurate to wi thin an<-absolute error of e = 1 G - 1 2 . S o l u t i o n : For convenience, let ENIU denote the integral error term EN,U (6 + 13i) in this example. The first task is to translate the absolute 31 Chapter 2. A Primer on the Gamma Function bound e = 1CT 1 2 on the error in computing T(7 + 13i) into a bound on EN^n. Let log(GPjv) denote the left hand side of (2.11), where G represents the true value of T(7 + 13i) and PN the product term. Let l o g G V . n denote the right hand side of (2.11) without the integral error term. Then (2.11) becomes \og(GPN) = \ogGN,n + E, N,n so that and the requirement is G = PN GN,n_eENtn _ GN,n N That is, it is desired that Pi N < e - 1 < from which \E N,n\ 2! 3! GN,U/PN\ < \GN,U/ PI N The infinite series is at most e in modulus provided Ex,n < 1, so that if \EN,U\ < e/(e|Gjv i T J /P/v|), the prescribed accuracy is guaranteed. This last estimate is self-referencing, but GN,U/PN can be estimated using Stirling's formula (2.13). That is, we require I-EV.nl < ee 2 - 1 2 l T Z z + 1 / 2 5 x 10 -12 upon setting z = 6 + 13i. Now that a target bound for \EN,U\ has been determined, n and N must now be selected to meet the target. In Table 2.1 are listed bounds on |i?jv,n| as given by (2.15) for z = 6 + 13i and various combinations o f ( i V > ) . 32 Chapter 2. A Primer on the Gamma Function n N = 0 N = 1 N = 2 N = 3 N = 4 1 1.9 x l f r e 1.6 x 10" 6 1.3 x 1CT 6 1.1 x 1 0 - 6 9.7 x 10" 7 2 3.7 x 10" 9 2.8 x 10~ 9 2.2 x K T 9 1.7 x 1CT 9 1.3 x 1CT 9 3 1.9 x 1 0 " 1 1 1.3 x 1 0 - 1 1 9.1 x 1 0 - 1 2 6.4 x i r r 1 2 4.4 x 1 0 " 1 2 4 1.9 x 1 ( T 1 3 1.2 x 1 0 " 1 3 7.3 x 1 0 - 1 4 4.6 x 1 0 - 1 4 2.9 x 1CT 1 4 5 2.9 x i c r 1 5 1.6 x 1 0 - 1 5 9.3 x 1(T 1 6 5.3 x 1 0 " 1 6 3.1 x 1 0 " 1 6 Table 2.1: Upper bound on |£JV ) T I(6 + 13i)\ in Stirl ing Series We can see that N = 0 wi th n = 4 is sufficient to produce the desired accuracy. Using these values wi th z = 6 + 13i, and leaving off the error term in (2.11) then gives lo g r(7 + 13z) « logG 0 ,4 = (6 + 13i + 1/2) log (6 + 13i) - (6 + 13i) + ^log27r 4 B-2j ^ 2 j ( 2 j - l ) ( 6 + 13i) 2i- 1 = -2.5778902638380984 + i 28.9938056395651838 . Taking exponentials yields the final approximation r(7+ 13i) « G0,4 = -0.0571140842611710-2 0.0500395762571980 . Using the symbolic computation package Maple wi th forty digit preci-sion, the absolute error in the estimate is found to be |T(7 + 13i) — £0,41 : 2.5 x 10-15, which is w e l l within the prescribed bound. Observe from Table 2.1 that N = 4 wi th 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 iV were selected in this example. The question of how best to choose n and iV given z and absolute error e is not an easy one, and wi l l 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 wi l l 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 wi th (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 UN,n(z) ( ^ M cos (9/2)) (2n + 2)(2n + l)(z + N)2n+1 ' which is the right hand side of (2.15). Recall that 9 = arg (z + N). M i n i m u m o f \UN,n(6 + 1 3 i ) | 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 |C/jv,n| = |^ 7/v,n(6 + 13i)|. To find the minimum value of |f^v,n|, wi th iV fixed, it is sufficient to determine the n value at which the sequence of |c7/v,n| changes from decreasing to increasing, which cor-responds to the first n such that \UN,n/UN,n+i\ < 1- Using the explicit form of the Bernoull i numbers (2.14), the task then is to find the least n such that (2TT) 2(Z + AT)2 cos 2 (9/2) C(2n + 2) (2n + 2)(2n + l ) C(2n + 4) That is, we want the least n satisfying < 1 |7r(z + N ) cos (0/2)] < n* ' 3 1 1 + + C(2n + 4) 2n 2n2 J ((2n + 2) ' The square root term is at most when n = 1 and decreases to 1 rapidly, so that if n is chosen \ir\z + N\ cos (9/2)1 (2-16) n 34 Chapter 2. A Primer on the Gamma Function a minimal error should result. In the case of z = 6 + 13i, cos (0/2) = \Jl/2 + S/V205 and N = 0, (2.16) predicts a minimum error for n = 38. Using Maple the cor-responding error is found to be j £7o,381 < 1.3 x 10~ 3 4 . This bound is impressive, but bear in mind that the corresponding Bernoull i number required to achieve this is \B76\ « 8 x 10 5 0 . This is the best possible accuracy for z = 6 + 13z wi th N = 0; greater accuracy is possible as N is increased from zero. For example, wi th N — 4, (2.16) predicts a minimum error for n — 47, at which |^4,471 < 6.7 x 1 0 - 4 2 . Uniform Bound of | Ujy<n(z) \ A s 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 Stirl ing series. It is worthwhile to determine this l imit ing accuracy in the form of the best possible uniform error bound as a function of N > 0. B y 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^(Z)\ wi l l be worst possible where \cos(0/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 |E/o,n| g r ° w s without bound. Assume then that N > 1. Wri t ing cos (0/2) = y/(l + cos0)/2 and cosO = u/\/v? + v2 puts (2.5.2) in the form I cos (0/2)\2n+2\u + iv\2n+1 ^ ^ r i ( u 2 + v2)n/2 (ytf + v* + Uy+1 which is clearly minimized when u and v are minimized, that is, at v = 0 and u = N. In terms of z and 0, this says that \UN,TI(Z)\ is worst possible at z = 0 = 0. B y (2.16), n should then be chosen n « \TTN], 35 Chapter 2. A Primer on the Gamma Function and we have 2n+2 (2n + 2)(2n+l)N2n+l 2(2n + 2)!C(2n + 2) (27r)2"+2(2n + 2) (2n + 1) N2n+l from (2.14) 2(2n)! ( 2 7 r )2n+2_/V2n+l since ((2n + 2) « 1 2 y ^ m (2n) 2 " e~ 2 n ( 2 7 r ) 2 r i + 2 A ? ' 2 n + 1 from (2.13) -2TTJV upon setting n ~ 7 r 7 V . So selecting n ~ \TTN~\ results in a uniform bound which decreases rapidly wi th N. The problem remains, however, that the larger N, and hence n becomes, so too do the Bernoulli numbers required in the Stir l ing series. In practice, unless one is dealing wi th 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 c o + 2^ TTT + e(z) T{z + l) = (z + a)z+1/2e-{z+a)(2ir)1/2 k=l z + k (2.17) which is valid for Re(z + a) > 0. The parameter a is real, N = \a] — 1, c 0(a) = 1 and ck{a) is the residue of T(z + l)(z + a)-^+1^ez+a(27T)-1/2 at z = —k. Explici t ly , for 1 < k < N, this is <*(«) = 4 = f e ^ , ( - * + « ) f c - 1 / 2 e - f c + a • (2.18) 2n(k-l)\ 36 Chapter 2. A Primer on the Gamma Function Spouge's formula has the very simple relative error bound eg (a, z) e(z) es{a,z) T(z + l)(z + a)-( z + 1 /2)e^(2vr)- 1 /2 < (27r)a+V2 R e ( z + a) ' provided a > 3. Thus for z in the right half plane Re(z) > 0, es(a,z) has the uniform bound \es(a,z)\ < V^(27r) a+1/2 • (2.19) 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 r(7 + 13i) from Example 2.1: E x a m p l e 2.2. Estimate T(7+13i) accurate to wi thin an absolute error of e = 10~ 1 2 using Spouge's method. S o l u t i o n : A n absolute error bound of e < 10~ 1 2 means the relative error must be bounded by \es(a, z)\ < r(7 + 13i) 2Q-12g6+13i 2TT(6 + 1 3 z ) 6 + 1 3 i + V 2 by Stirling's formula = 1.3 x 10 - i i B y plotting |es(a, 6 + 13z)|, 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 Cfc(a) 0 +1.0000000000000 x 10° 1 +1.3355050294248 x 10 5 2 -4.9293093529936 x 10 5 3 +7.4128747369761 x 10 5 4 -5.8509737760400 x 10 5 5 +2.6042527033039 x 10 5 6 -6.5413353396114 x 10 4 7 +8.8014596350842 x 10 3 8 -5.6480502412898 x 10 2 9 +1.3803798339181 x 10 1 10 -8.0781761698951 x 10" 2 11 +3.4797414457425 x 10~ 5 12 -5.6892712275042 x 1 0 " 1 2 Table 2.2: Coefficients of Spouge's Series, a = 12.5 then gives T ( 7 + 13i) « -0.0571140842611682 - i0.0500395762571984, which differs in absolute value from T(7 + 13i) by less than 2.9 x 10~ 1 6 . 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 TV is straightforward. • Spouge's work is important for several reasons. The first is that the coefficients Cfc(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. A n d finally, Spouge's approximation and error estimates apply not only to T(z + 1), but also to the digamma function ty(z + 1) = d/dz [logr(z + 1)] arid tr igamma function ty'(z). To see the link between (2.17) and Lanczos' formula (1.1), write a = r + 1/2 and resolve the first terms of the series (1.2) into partial 38 Chapter 2. A Primer on the Gamma Function fractions - a 0 ( r ) + a i ( r ) z + 1 z H \-aN(r) z~-(z-N + l) (z + l)---(z + N) + e{z) N h(r) z + k Comparing this wi th (2.17), the bk(r) obtained from truncating the Lanczos series are the approximate residues of T ( 2 ; + l ) ( 2 ; H - a ) _ ^ 2 + 1 / ' 2 ^ e 2 + a ( 2 at z = — k, and the larger N becomes the better the approximation. The standard definition and properties of the gamma function can be found in many works dealing wi th 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 A r t i n [2] is of interest for the treatment of T(x) for x veal 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. . 2.7 Additional Remarks 39 Chapter 3 The Lanczos Formula We arrive at last at the examination of Lanczos' paper itself, beginning in this chapter wi th the derivation of the main formula (1.1). The derivation consists of three steps, which in a nutshell can be broken down as /•CO r(z + l / 2 ) = / t'-^e^dt Jo x ^ [v(l-logv)}z-1/2vrdv Jo (step I) z + r + l/2)z+1^e^z+r+1^V2 V2v{6)rsm6 r/2 x / co s 2 2 0 J-n/2 \ogv(6) dO (step II) = {z + r + l/2)z^2e^z+r+l^V2 / cos 2 2 0 J —7T X T / 2 ^ + f > ( r ) c o s ( 2 A ; 0 ) fe=0 dO . (step III) 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 wi th 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 wi l l be standard throughout the re-mainder of this work: D e f i n i t i o n 3 .1 . Define H0(z) = 1 , and for k > 1, T(z + l)T(z + l) k K ) T{z-k + l)T{z + k + \) 1 ~ (z + l)k(z + l)-k _ z • • • (z - k + 1) ~ (z + 1) v (z + fc) ' Using this new notation, (1.1) may be restated 1: 0 0 • r ( z + 1) = (z + r + l / 2 ) 2 + 1 / 2 e-(z+r+1^ ^ ak(r)Hk(z) . (3.1) fc=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 = o + it where o = Re(z) and t = Im(z); 2. For a, b G . C wi th a not a negative real number or zero, define ab = e 6 1 o s a where Im(loga) is between — TT and TT. 3.1.1 Preliminary Transformations Beginning wi th pOO T(z + 1)= fe-'dt, Re(z) > - 1 , (3.2) Jo make the substitution t —> at where R e ( a ) > 0 to obtain pOO T(z + 1) = a z + 1 / tze~atdt . Jo The validity of this substitution for complex a can be shown using contour arguments or via the following lemma 2 : L e m m a 3.1 . For Re(a) > 0 and Re(z) > —1, poo T{z + 1) = a z + l / tze~atdt . Jo P r o o f o f L e m m a 3.1: For fixed z wi th Re(z) > —1 and a > 0 real, replacing t wi th at in (3.2) gives poo T{z + 1) = a z + l / tze~atdt . Jo Viewed as a function of a, the right hand side is analytic in the right half a-plane, and equals the constant T(z + 1) along the positive real 2 This result can be found as an exercise in [30, p.243]. For real a > 0 the result is elementary, however Lanczos makes no comment as to its validity for complex a. 42 Chapter 3. The Lanczos Formula axis. Hence the integral must equal T(z + 1) throughout the right half a-plane. • Next set a = z + r + 1, where r > 0, so that roo T(z + 1) = (z + r + l ) z + 1 / tze~(z+r+1)t dt, Re{z) > - 1 . (3.3) Jo Here r > 0 is a free parameter which plays a key role later. Now reduce the bounds of integration to a finite case v ia the change of variable t = 1 — logTj, dt = (—l/v)dv: r(z + l ) = (z + r + l)z+1 / (l-\ogv)ze-{z+r+1){1-]ogv)—dv Je v = (z + r + i ) * + i e - ( * + H - i ) f (i _ \ogv)zvzvrv-dv Jo v = (z + r + i ) ^ + i e - ( 2 + r + 1 ) / [v(l - log v)]z vrdv . (3.4) Jo The integrand in this last expression is 0(va+r~e) for every e > .0 as v —> 0, while it is 0(\v — e\a) as v —> e . Equation (3.4) is therefore valid for Re(z) > — l . 3 The substitution t — 1 — log t> appears to be the key to the Lanczos method in the sense that it puts al l the right pieces in all the right places for the series expansion to come. It has the effect of peeling off 3There is an error in [17, p.30, Eq.(l)] at this point concerning the domain of convergence of the formula. There the author makes the replacement z —> z — 1/2 and states (using a to signify what is here denoted r): F(z + 1/2) = (z + a + l / 2 ) 2 + 1 / 2 e x p [-(* + a + l /2)]F(z) , F(z) = JQ[V{1 -\ogv)]z-1/2vadv, Re(z + <r + 1/2) > 0 . However, the integral F(z) diverges with, for example, z = —3/4 and a choice of o = 5/4. To see this, note that v(l — logv) < (e — v) on [0, e], so that [v/(e - v ) ] 5 / 4 < K 1 - log v ) ] _ 5 / 4 u 5 / 4 , whence [e[v/(e - v)f4 dv < fe[v(l - l o g 7 j ) ] - 5 / 4 7 j 5 / 4 dv . Jo Jo But the left hand side /Q e[v/(e — ? J ) ] 5 / 4 dv = oo. 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 a —» (z + r + 1) and t —> 1 — logf , 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{6) 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(l — logv) is shown in Figure 3.1 so that v may be implici t ly defined as a function of 9 v ia cos 2 9 = v(l — logv) , where 9 = —ir/2 corresponds to v — 0, 8 = 0 to v = 1, and 9 = ir/2 to v = e. 0 1 v Figure 3.1: Pa th of v(l — logv) Making the substitution v(l—logv) = cos2 9, dv = 2 sin0 c o s 0 / l o g f d9 44 Chapter 3. The Lanczos Formula in (3.4).then gives r'z +1) - (z + r + i)-+i e-(-+-+D f / 2 (cos 2 ey vr 2 s™9 0 0 8 0 d 9 J-n/2 loEV = (z + r + i)-+i e-(-+-+D f / 2 C o s 2 z + 1 0 2 ^ 8 i n ^ dO . (3.5) J-ir/2 l(>gV The integrand in (3.5) is 0 ( | 0 + 7 r / 2 | 2 f f + 2 r + 1 ) near 0 = -TT/2, 0 (1) near 0 = 0, but is O( |0 — 7r/2| 2 c r + 1) as 0 —> 7r/2 so this expression in once again valid for Re(z) > — 1. Moving a A / 2 outside the integral and replacing z wi th z — 1/2, the reason for which wi l l become apparent later 4 , yields / .TT /2 T(z +1/2) = Pr(z) / c o s 2 z 0 ' - T T / 2 where Pr(z) is defined Pr(z) = V2(z + r + i / 2 ) 2 + V 2 e - ( 2 + r + i / 2 ) _ Now denote by / r (0 ) the term in square brackets in (3.6) and by fE,r{9) its even part [fr(9) + fr{—9)] / 2 . Noting that the odd part of the inte-grand integrates out to zero finally puts (3.6) in the form rn/2 T{z +1/2) = Pr(z) / cos2z 9 fEAO) d9 (3.7) for Re(z) > —1/2 and r > 0. Equation (3.7) is the starting point for the series expansion in (3.1). 4Shifting the argument to 2 — 1/2 at the outset, that is, starting the derivation with roo r ( z + l / 2 ) = / t'-^e^dt, R e ( z ) > - l / 2 Jo 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. y/2vr sin 0 logv d9 , (3.6) 45 Chapter 3. The Lanczos Formula 3.1.3 Series Expansion of fE,r The next step in the development is the replacement of JE,r(9) by a series which can be integrated term by term. Lanczos first constructs a Taylor series in sin# for fr(0) 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 fs,r w i l l be examined later, but for the time being assume that this function may be represented by a uniformly convergent Fourier series on [—TT/2, 7r/2]. Further, this series wi l l contain only cosine terms since fs,r is even. Thus fE,riP) = 9 a o ( r ) + ^2ak(r) C O S 1 fc=l 0 0 = ^jT^' otfc(r) cos (2k0) (3.8) fc=0 where the coefficients ak(r) are given by *k(r) = - fE,M c o s (2fc0) d9 . n J-n/2 O n 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 T(z + 1). This identity is stated as a lemma, the. proof of which can be found in [17, p.16]: L e m m a 3.2. For Re(z) > - 1 / 2 cos 2 2 9 cos (2k9) dO = ^ 1 . ' Hk(z) (3.9) TT /2 r ( z + 1) Picking up where we left off wi th (3.7), and using (3.8) to replace 46 Chapter 3. The Lanczos Formula fE,r(Q), we obtain T(z +1/2) = Pr(z) cos2z9 ak(r) cos (2k6)d8 ' - T T / 2 fc=0 = Pr{z)^2 ak{r) I cos2zdcos{2k9)dd k=0 "~7r/2 oo = fi(^)v^r^ + 1(f £ ' a f c ( r ) f f f c ( * ) . (3.10) ^ ^ fc=0 Now it is just a matter of eliminating r(z + 1/2) on both sides of (3.10), the reason for replacing z wi th z — 1/2 in (3.5), and moving the T(z + 1) to the left hand side to get the final form: oo Y{z + 1) = (z + r + l/2)z+^2e-^+r+1^V2?} £ ' a k ( r ) H k ( z ) , (3.11) k=0 where again, Re(z) > —1/2, and r > 0. We wi l l 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 wi th equation (3.4), r(z + l ) = (z + r + l ) z + 1 e - { z + r + 1 ) f [v{l -logv)]zvrdv , Jo this time define v(x) implici t ly with 1 -x2 = v(l -logv) , (3.12) where x = — 1 corresponds to v = 0, x — 0 to v = 1, and x = 1 to v = e. This yields T(z + 1) = (z + r + i ) * + i e - ^ + r + 1 ) j1 (1 - x2f '2xvr logv dx . (3.13) 47 Chapter 3. The Lanczos Formula Once again let Pr(z)= V2 (z + r + l/2)z+ll2 exp [-{z + r + 1/2)], re-place z wi th z—1/2, and denote by JE,r{x) the even part of \[2xvr j logt>, giving f1 fir T(z + 1/2) = Pr(z) J i l - x2)z fE,{x) - j = = . (3.14) The motivation this time for the replacement z —> z—1/2 is a bit clearer: wi th the introduction of the weight function l / \ / l — x2, the integral i n (3.14) reveals itself as an inner product on £ 2 _ 1 ^(dx/y/l — x2) which has the Chebyshev polynomials as orthogonal basis. Expanding fE,r{x) in a series of even Chebyshev polynomials (the justification for which wi l l follow) produces 0 0 / E , r ( a : ) = E 'c f c (r)T 2 f c (x) (3.15) r1 T(z + 1/2) = Pr(z) / (1 - x2)z ck(r)T2k(x) 1 i—n X fc=0 so that dx Ck{r)T2k{x) — fc=0 which upon integrating term by term gives T(z +1/2) k=0 vT 2 ' OO r (z+ 1/2) = PM^r^l™ ^ Ck{T){-l)k Hk{z) . The term by term integration uses an equivalent form of the iden-•tity (3.9): . Finally, mult iplying through by r(z + l)/r(z + 1/2) and writ ing ° fc ( r ) = (— 1) f c cfc( r") once again yields T(z + l) = (z + r + l/2)z+1/2e-^+r+1^V2^Y^ak(r)Hk^) • (3-16) fc=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 T(z+l/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 fs,r be equal to a uniformly convergent Fourier (resp. Chebyshev) series is unnecessary, if we ask only that fE>r be square summable wi th respect to Lebesgue measure. Take for example the Fourier setting. Parseval's theorem [24, p.91] says that for even / , g £ L2[—7r/2,7r/2] wi th g real, 2 r / 2 w . where fk denotes the Fourier coefficient fk = - f{9) cos (2k9)d9 . In the present case, the Fourier coefficients of p(9) 2 r / 2 Pk = - cos2z (9) cos (2k9)d9 T J-n/2 _ 2 r(z + l/2) r(z + i) H k { z h from Lemma 3.2, while those of fE,r are ak(r), assuming fEtT £ L2[—TT/2, n = cos 2 z(#) are 49 Chapter 3. The Lanczos Formula Thus by Parseval's theorem, (3.7) becomes simply / •TT /2 T(z + 1/2) = Pr(z) / cos2* 0 fE,r(0) d9 ' - T T / 2 oo 71" ^ — s, i = p r ( ^ 2 2 ^ Pkak(r) k=0 OO = P r i z ) ^ ^ ^ £ak(r)Hk(z) , 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 0 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-sibili ty 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 fo{0) 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 wi l l 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 . A r m y 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 wi l l 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 a^(r) (or equivalently Cfc(r) in the Chebyshev set-ting). There are several methods for computing the coefficients and these wi l l be examined in due course. For the moment, note Lanczos' own formulation: the constant term of the series is Here J>(j) = 2- 1/ 2r(j + 1/2)(j + r + l / 2 ) - ^ 1 / 2 exp (j + r + 1/2) and C2j,2k is the coefficient of a;2-7 in the 2kth 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 wi l l 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 afc(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). (3.17) while for k > 1 ^ C2j,2k^:r(j) • 3=0 (3.18) 51 Chapter 4 The Functions v and JE,T The derivation of Lanczos' main formula (1.1) in Chapter 3 relies on the validity of expanding the implici t ly defined function fE,r(Q) of equa-t ion (3.7) (resp. fE,r{x) °f (3-14)) into an infinite interpolating Fourier (resp. Chebyshev) series. Furthermore, the smoothness and summa-bil i ty properties of f E < r determine the asymptotic growth rate of the coefficients ak(r) (resp. Cfc(r)), and consequently the convergence of the series (1.2). In this chapter the properties of the implicit functions v, fr and f E > 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) e C''-r-'[—1,1] and is of bounded variation, so that the expansion (3.15) is justified. (ii) In the Fourier setting, fE,r(9) G C'- 2 r J [—7r/2, TT/2] and is of bounded variation, thus justifying the expansion (3.8). (iii) ff} 6 L 1 [ - T T /2 , T T /2] , where n = \2r]. (iv) ak(r) = 0 ( A r r 2 r l ) as -> oo. 4.1 Closed Form Formulas We begin by finding explicit expressions for v and fE,r in terms of Lambert W functions Wo and W-\. 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 jE%r smoothness properties of v and fE,r can be deduced directly from those of Wo and Th i rd , the expression for fEt1. 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) o n [—1> !]• These formulas are in terms of the two real branches W_\ and Wo of the Lambert W function defined implici t ly by W(t)ew® = 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 implici t ly v ia 1 ~x2 = v{l -logv) , (4.1) where x — — 1 corresponds to v = 0, x = 0 to v = 1, and x = 1 to v = e. Lett ing w = log (v/e), this is equivalent to * 2 ~ 1 = we whence so that V - 1 w = W v fx2 - I - — exp W e w(^) (4.2) For — 1 < ^ < 0 , 0 < W < 1 which corresponds to the real branch W-i of the Lambert W function. For 0 < a : < l , 0 < t ; < e which corresponds to the branch Wo. Using the Heaviside function H(x), and letting y(x) = [x2 — l ) / e , v can thus be expressed as the single formula vM = H ( _ x ) y(x) e W-Mx))H{ X ) + Wo(y(x))H{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-\ (y{x)) —> —oo, 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. \ i . X1 -1 JL ~~"~e a -+~ \ I j ] — t I | l j 1 ... L I I | i I 1 1 0  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) = \f2xvr/logv, note that logv = 1 + W((x2 l ) / e ) from (4.2), so that fr(x) may be expressed fr{x) = V2 X l + W-i (y(x)) H(-x) ( v{x) V \W0(y(x))J 1 + W0(y(x)) H(x) 54 Chapter 4. The Functions v and fE>r while its even part fE,r(x) = ifr(x) + fr(—x)] /2 becomes IEAX) _ N_ ( v(x) V \Wo(y(x))J 1 + W0(y(x)) l + W-xiyix)) W-i(y(x)) (4.4) Plots of V2fE,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) — y/2 xvr/ logv that y/2 / £ i r ( ± l ) / e r = 1 if r > 0 and is +oo 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-i I I i • i i i i | I i I_ | 1 ' r r= C ) • — j — i | | ! — | |... \ \ 1 i ! 1 ! i \ i i. > ! ! 1 I I 1 I I 1 f | r -= 1 / / ' r -_ r / | r -= I i i 1 i 1 i - 1 0 1 Figure 4.2: V2 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 /^,r(cc) The smoothness of vr, fr and fE>r can be seen to be essentially the same problem wi th the observation that y/2 fr(x) = d/dx[vr+l/(r +1)]. This follows from the fact that 1 — x2 = v(l — logv) from which d , . 2x , In this section it is shown that vr, fr (and hence fE,r) have |_rJ contin-uous derivatives on [—1,1]. Smoothness on (—1,1] Away from x = — 1 (where v > 0), the smoothness of ?/ is precisely the same as that of v. A s such, it is enough to consider v(x) on (—1,1]. O n ( - 1 , 0) U (0,1), the functions Wk((x2 - l ) / 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 = e x p W 0 ( ( ^ 2 — l ) / e ) which is again smooth since Wo(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. B y expressing the right hand side of (4.1) as a Taylor series, it is not difficult to show that (4.1) may be writ ten 1 1/2 x = (v — 1) 5( 1 > ( j + l ) ( 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 4>(v, x) = x-h(v). Now 0(1,0) = 0, d(f)/dv\(h0) = -l/y/2, and 0 € C°° at (1,0). Therefore by the implicit function theorem there are open sets U C R2 and W C M and a function g : W - • R such that g e C°°, (1,0) e U, 0 e W, g(0) = 1, and <f>(g(x),x) = 0 for every x e W (see [19, p. 103]). Thus v = g(x) is smooth about x = 0. Thus v is smooth as a function of a; on (—1,1], from which vr 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. A t this point the role of the parameter r becomes crucial and must be taken into account. W h e n r = 1, v and v' = 2x/\ogv are zero as x —» — 1 + , while v" = 2 / l o g v — 4x2/[v(logv)3] —> co there. Thus v is once continuously differentiable at x = —1. Raising v to the power r > 1, on the other hand, forces vr to zero more rapidly as x =—> — 1 + , and so has the effect of smoothing out the singularity there. Begin wi th a simple calculus result: if r > 0 and k < 0 then l im | x J V' ( log i ; ) f c | = l im vr |logi>| f c x—*—1+ v—>-0+ < oo . (4.6) W i t h this in mind, consider the smoothness properties of the gen-eralized function g(x) = xivr(logv)k where j and k are integers and r is real. Using (4.5) to differentiate g(x) yields [z-VQogtOkl = dx J jxj-1vr(logv)k + 2rxj+1vr-\logv)k-1 + 2kxj+1vr-\\ogv)k-2 which can be abstracted as G(j, r, k) V$ jG(j-l, r, k)+2rG(j+l, r - 1 , k~l)+2kG(j+l, r - 1 , k-2). (4.7) In the present case, vr corresponds to,G(0, r, 0), while fr corresponds to 2 G ( l , r , —1). Now beginning wi th G(0 , r , 0) and differentiating re-peatedly, (4.7) shows that wi th each derivative, the k component of any resulting G terms remains zero or negative while the least r com-ponent reduces by 1. B y (4.6), at least the first |j-J derivatives of vr are bounded as x =—>• — 1 + . The case r = 1 shows that this bound is best possible. Since fr(x) = d/dx[vr+1 /(r + 1)], it follows that fr has [rj derivatives as well. In summary 57 Chapter 4. The Functions v and fE,r Theorem 4.1. (i) fr(x) andvr are smooth on (—1,1]; (ii) fr(x), fE,r(x)> a n d v r are m C ^ [ - l , l ] . Before proceeding to the examination of fr as a function of 0, we end this section wi th 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]. P r o o f of Theorem 4.2: It is enough to show f'r{x) > 0 on (—1,1). Since logv and vr is non negative and increasing, it is enough to show that x h(x) = logv is non negative and increasing. It is clear that h(x) = v'(x)/2 > 0. Now , logv - 2a; 2/(vlogv) h (x) = j- , ( logv) 2 so it is enough to show that the numerator p(x) = logv — 2x2/(v logv) _ v(log (v))2 - 2 + 2 v - 2 v l o g (y) 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(l) — / r (—1), from which 58 Chapter 4. The Functions v and fs,r C o r o l l a r y 4 .1 . fr(x) and fE,r{x) a r e functions of bounded variation on [-1,1]. From the continuity of fE,r{x) a n d Corollary 4.1 it follows that fE,r{x) m a y be expressed as a uniformly convergent Chebyshev series on [-1,1]. 4.3 Smoothness of v r(0), /P(0) and fE,r(e) B y setting x — sin 9 and d dx d d~9 = d9dx ' it follows easily from the previous section that vr(9), fr(9) and fE,r(9) have at least [rj continuous derivatives on [—TT/2,TT/2]. In fact, this argument shows that these functions are smooth on (—TT/2, TT/2], wi th a singularity at 9 = — TT/2. In this section, the precise order of the singularity of fr(9) (and consequently fE,r(9)) at 9 = —7r/2 is determined. This wi l l 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 G e n e r a l i z e d fr(0) Define G(/3,j, k, I; 9) as G((3, j, k,l;9) = v>3 sinj 9 cosk 0(log v)1 , (4.8) where 8 e R , j, fc, I E Z . Note that fr(9) = V2G(r, 1,0, - 1 ; 9). From tedious application of derivative rules and using the fact that 2 sin 9 cos 9 59 Chapter 4. The Functions v and fE,r G((3, j, k, I; 6) has the derivative form d de [G(f3,j,k,l;0)} = 2 / ? G ( / ? - U + l,fc + l , Z - l ; 0 ) +jG(f3,j-l,k+l,l;0) -kG((3,j + l , k - 1,1; 6) +2lG((3-l,j + l,k + l,l-2;0) . (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 wi l l help to clarify the dependence of G on (3, j, k and /. Begin wi th an easy result from calculus: Lemma 4.1. Suppose (3 G R and I G Z . Then l i m v13 IlogvI ' 0 if (5 > 0 oo ifp<0 = < 1 «/)9 = 0 , / = 0 oo if(3 = 0, I >0 0 i / 0 = O, Z < 0 . From this follows Lemma 4.2. Suppose (3 G R and j , k,l £Z. Then l i m G((3,j,k,l;6) = { 2 0 if(3 + k/2>0 oo if(3 + k/2<0 (-l)l+j if (3 +k/2 = 0 and k/2 + l = 0 ( - l ) ^ ' - o o if(3 + k/2 = 0 and k/2 + l > 0 0 if(3 + k/2 = 0 andk/2 + l <0 Proof of Lemma 4.2: Near 0 = -ir/2 cos0 = ? j 1 / 2 | l - l o g v | 1 / 2 Also recall that as # —> — n/2 , v —> 0 + . Now write 60 Chapter 4. The Functions v and fs,r l im G(8,j,k,l;9) l im v^sin- 7 #cosfc#(logt>) 6 e = l im v0vk/2\l-\ogv\k/2smjO\\ogv\l(-l)1 = l im ^ + f c / 2 | l o g v | ' + f c / 2 ( - l ) ^ 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 —n/2. To apply the results of the previous section to fr, abstract the process of differentiating G(/3,j, k, I; 9) as follows. Identify wi th G(8,j, k, I; 9) the monomial x/3yjzkwl, and view differentiation of G(8,j, k, I; 9) as the linear mapping : Higher order derivatives of G(8,j,k,l;9) are then obtained recur-sively as D2 = D o D, D3 = D o D o D, etc. For consistency, let D° denote the identity mapping. 4.3.3 Application to fr We arrive finally at the analysis of fr. Theorem 4.3. fr(0) is continuous, increasing and hence of bounded variation on [—TT/2, TT/2}. 4.3.2 Derivatives of G(f3, j , fc, I; 6) D : xl3yjzkwl i-> 28xP-lyj+1zk+1w1-1 +jx0yi~1zh+1wl -kx^yj+lzk-lwl +2lxP-lyi+lzk+lw1-2 . (4-10) 61 Chapter 4. The Functions v and fE,r P r o o f of Theorem 4.3: This follows from Theorems 4.2 and 4.1 upon writ ing x = sin 9 and dfr dfr dx lO = dxdO ' Theorem 4 .4 . fr(9) has [2r\ continuous derivatives at 0 = —TT/2. P r o o f of Theorem 4.4: Using the notation of Section (4.3.1), write fr(9)/\/2 = G(r, 1,0, —1;0). B y Lemma 4.2, fr is continuous at 9 = —TT/2 for r > 0. Now consider the associated sequence D°(xry1z°w~1), Dl{xrylzQw-1), D\xrylzQw-1), . . . (4.11) The task is to determine the least n such that Dn+1(xry1z°w~1) contains a term whose powers meet one of the oo conditions in Lemma 4.2. G(f3,j,k,l;6) w i l l then be n-times continuously differentiable at - 7 T / 2 . Observe first that the third and fourth conditions of Lemma 4.2 (P + k/2 = 0 and k/2 + l > 0) can never be satisfied by a term of one of the Dn(xry1z°w~1) in the sequence (4.11), since these conditions imply j3 — I < 0, which translates into deg (x) — deg (w) < 0 in some term of Dn(xry1z°w-1). Bu t deg (a;) - deg (w) = r + 1 > 0 i n D ° and (4.10) implies this quantity is non-decreasing wi th n. The only other oo condition in Lemma 4.2 is satisfied when deg (x) + deg(z) /2 < 0 for a term in some Dn(xry1z°w~1). In D ° , deg(x) + deg(z) /2 = r > 0. Equation (4.10) implies that the least value of deg (a;) + deg(z) /2 taken over the summands of Dn(xry1zQw~l) decreases by at most 1/2 wi th each increase in n . Further, since deg(w) = —1 in D°, the fourth term of (4.10) guarantees that this decrease wi l l indeed occur wi th each n. Thus, it is enough to determine the largest integer n such that r — n(l/2) > 0, which is the largest n such that n < 2r, i.e. n = [2r\. In other words, for r > 0, fr(9) has |_2rJ continuous derivatives. Lemma 4.2 also gives that al l of these |_2rJ derivatives are zero. Theorem 4 .5 . /rV/2) = 0 for n a positive odd integer. 62 Chapter 4. The Functions v and JE,T P r o o f of Theorem 4.5: Aga in identify fr(9)/y/2 wi th G(r, 1, 0, - 1 ; 9) and consider the associated sequence of derivatives D°{xrylz°w-1), D1(xry1z°w-1), D2(xrylzQw~l), ... B y (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 cos0 which is zero at 9 = TT/2. Theorem 4 .6 . / r ( n ) e L1[-TT/2,TT/2] where n = \2r]. P r o o f of Theorem 4.6: It has been established that / is smooth on (—TT/2, TT/2], S O it is sufficient to bound the growth of/r"^ as 9 —» —TT/2+. A g a i n it is easiest to work wi th the generalized fr as given in Section 4.3.1, from which it is not difficult to determine conditions on 8,j,k and I which guarantee summability of the integral. Beginning wi th / = /"^ \G(8,j,k,l;9)\ d9 T / 2 / .TT /2 • / - T T / 2 v13 sin-7 9 cosk 9 (\ogv)1 d9 change to v variables wi th the substitution v(l — logv) — cos 2 9: v " [ l - v(l - logu)p'/2[u(i - \ogv)]k/2 logv dv j0 | ( logu) ' [ l - V(l - l0gu)] 1/2[ v(l - logv)]1/2 - v{l - l ogv ) ]^ - 1 ) / 2 [v ( l - logv)f-^2 \logvl1-1 Jo dv Asv^ 0 + , the integrand is O^+^-^^ogv)^2-^1/2). Provided 8 + k/2 - 1/2 > - 1 , that is, 8 + k/2 + 1/2 > 0, the integral converges. In the case at hand of fr(0) — \/2G(r, 1,0, — 1;9), the exponent 8 + k/2 — 1/2 starts off at r — 1/2 and decreases by 1/2 wi th each derivative as (4.10) shows. Thus if n is chosen to be the largest integer such that r - l / 2 - n ( l / 2 ) > - 1 63 Chapter 4. The Functions v and fE<r then convergence of J^^lf^ldO is guaranteed. That is, if n is the largest integer such that n < 2r + 1, which is just [2r], then / r " ' G L 1 [ -T T /2 , TT/2]. • The results of Theorems 4.3, 4.4, 4.5 and 4.6 extend easily to the even part of fr: C o r o l l a r y 4 .2 . Suppose r > 0 and f , m fr(0) + fr(-0) jE,r\P) = 2 ' Then (i) fE,r is a continuous function of bounded variation on [—7r/2,7r/2] such that fE,r(—7r/2) = fE,r(^/2) = e r / v / 2 and fE,r(0) = 1; (U) fE,r e Cw[-7r/2,7r/2], w/ieren = L2rJ, a n d ( - T T / 2 ) = / g ( 7 r / 2 ) 0 / o r n < [2r J an odd integer; (tit) fEn} G L 1 [ -T T /2 , TT/2], w/iere n = \2r]. 4.4 Fourier Series and Growth Order of Coefficients We conclude this chapter wi th the application of Corol lary 4.2 to the justification of the replacement of fE,r by its Fourier series i n equa-t ion (3.8), and to the estimate of the growth order of the Fourier coef-ficients ak(r) i n equation (1.1). First of al l , from (i) of Corollary 4.2 it follows that fE<r has a Fourier series expansion which converges uniformly on [—n/2, TT/2] (see [25]). Since fE,r(9) is even, this series wi l l contain only cosine terms. Thus ^ oo 1 fc=l 64 Chapter 4. The Functions v and JE,T for coefficients *k(r) = - f ' 2 fEAO)cos(2k0)d8 * J-n/2 2 r ' 2 T / 2 B y substituting 9 = TT/2 and 9 = 0 at this stage we get er 1 0 0 4T? = ^ o ( r ) + £ ( - l ) f e a f c ( r ) and 1 0 0 l = ^ o ( r ) + £ a f c ( r ) , (4.12) useful facts for testing the accuracy of computed afc(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. A s k -> oo, a f c(r) = 0(k~^). P r o o f of Theorem 4.7: Beginning wi th 2 r / 2  fn/ ak(r) = - / fEAO) cos (2k9)d9 K J-n/2 integrate by parts n = \2r \ times to get 2 , ( » ) Here G f e ( r ) = / T 2 f{^{9) {cos (2^} 1 sin {2ke)}n de • ( 4 i 3 ) r / o m M • / o / m i / cos (2A;6>) if [2r"| is even, {cos (2*0) | sm (2k9)}n = j ^ ( ^ . f ^ . g ^ From (4.13) it follows that 1 2 r / 2 T / 2 whence ak{r) = 0(k-W) . (4.14) d9 , 65 Chapter 5 Analytic Continuation of Main Formula A t the end of the derivation of the main formula (1.1) i n Section 3.1, it 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 F(z + 1) in this region. Begin wi th some notation: Def in i t ion 5.1. For r > 0, define the open set nr = {zeC\ Re(z) > -r and z ^ - 1 , - 2 , - 3 , . . . } . 5.1 Domain of Convergence of Sr(z) Theorem 5.1. Suppose r > 0. Then the series oo Sr(z)= J^ak(r)Hk(z) converges absolutely and uniformly on every compact K C f L . Conse-quently, the series defines an analytic function on Qr. P r o o f of Theorem 5.1: First, observe that each term of the series is analytic on i}r. Now suppose K is a compact subset of i1r, let D(0, R) 66 Chapter 5. Analytic Continuation of Main Formula be a closed disk of radius R containing K, and let S be the distance between K and C \ fir. It wi l l first be shown that \Hk(z)\ = O (k2r~2S~l) on K. For k > R, by the reflection formula (2.10), T(z + l)T(z + l) Hk(Z) = T(z + 1 + k)T(z + 1 - k) (-l)k+1[T(z +I)}2 sin (nz) r(l-z + fc) Tr r(l+ * + *;)(*;-*) , r(i-z + fc) ^ru+ * + *)(*-*) ' say-Now h(z) is analytic on i1r, so 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, r(i - z + k) T(l + z + k)(k - z) 2nez-k(k- z)k-"+1^(l + Si(k)) 1 2ire-z-k(k + z)k+*+V2{l + 62(k)) (k - z) „22 k — z~ k + z ( k - z \ 1 / 2 1 (l + fr(fc)) {k-z)(k + z)\ \k + z) (k - z) (1 + 62(k)) where Sj(k) —> 0 as k —» oo and the 5j do not depend on z. Selecting 67 Chapter 5. Analytic Continuation of Main Formula z G 73(0, R) to make this last expression as large as possible, r(i -z + k) T(l + z + k)(k- z) < e 2R 'k + R \ K fk + R \ 1 / 2 ( l + frflfc)) 1 k-RJ \ k - R j (1 + 62(k)) (k - R)*°+i < e 2R 'k + R \ K fk + R \ 1 / 2 (l + Si'h)) 1 k-RJ \ k - R j (1 + S2(k)) (k - R)2(-r+s)+1 = o(k2r~2S-1) . This last inequality follows since, as k —> oo, „2R k + R \k fk + R \ 1 / 2 (l + S.jk)) 4 R k-RJ [k-RJ (l + 82(k)) ~"e Recalling that \h(z)\ is bounded, the conclusion is \Hk{z)\ = O (k*-**-1) on K. Since |a f c(r) | = 0(k~^) from Theorem 4.7, it then follows that \ak(r)Hk(z)\ = 0 ( A > - r a r l - i - M ) = o ( r 1 - 2 < 5 ) , whence, for some positive constant A and any n > 1 £ K ( r ) i 7 f c ( z ) | < A £ — ^ < o o , completing the proof. k=n 68 Chapter 5. Analytic Continuation of Main Formula 5.2 Analytic Continuation Theorem 5.2. For r > 0, the function oo V2^{z + r + l/2)z+1/2 e-{z+r+l'2) ak{r)Hk{z) is analytic on ilr. P r o o f of Theorem 5.2: The factor V2n (z + r + l / 2 ) 2 + 1 / 2 is analytic since Re(z + r + l / 2 ) > 0 on O r . The factor exp [-(z + r + 1/2)] is entire. The series factor is analytic by Theorem 5.1. Theorem 5.3. oo T{z + 1) = \Pht [z + r + l / 2 ) 2 + 1 / 2 e - ( z + r + l ' 2 ) E' ak{r)Hk{z) k=0 on Clr. P r o o f o f Theorem 5.3: The right hand side is analytic on Q r and equals T(z + 1) on Re(z) > —1/2. Hence by the principle of analytic continuation, it equals T(z + 1) throughout i}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 wi th the bound given following equation (3.11). However, the proof of Lanczos' bound eludes me. Observe that Sr(z) may be expressed Sr(z) = T{z + 1 ) ( 2 T T ) - 1 / 2 ( Z + r + i / 2 ) - ( ^ 1 / 2 ) ^ + r + i / 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 on the half plane Re(s) > C\, and f(s) is expressible as a series f(s) = J2T=i Qkk~a on Re(s) > c2 > Ci, where the coefficients qk are eventually positive. Then the series actually converges on the larger set Re(s) > c\. The best that can be achieved wi th the techniques used here is convergence on Re(z) > - |~2r] /2 , a slight improvement over Theo-rem 5.3. If the bound a^(r) = 0(k~^) could be improved to a^(r) = 0(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)k2r+l\ is compared against — log \ak{r)k^ for the first 150 coefficients wi th r = 1.25 (— log |ajt(r)| is also plotted for reference). O f the two lower curves, the flatter one corresponding to — log \ak(r)k2r+l\ indicates the tighter bound. The difficulty lies in finding a more precise bound on the Fourier coefficients than the rather cumbersome estimate ak{r) = 0(k~^). To do this requires a better understanding of the asymptotic behaviour of fr(0) about the singularity 0 = —TT/2. Lett ing h = (0 + n/2)2, it is true that as 0 —> —7r/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 Chapter 6 Computing the Coefficients In order to use the formula (1.1) in practice, the Fourier coefficients afc(r) must first be determined. The method given by Fourier theory, namely the direct integration i n not practical due to the complicated nature of the factor fE,r(0) in the integrand. There are, however, a number of other methods for computing the afc(r), and in this chapter these various methods are discussed along wi th 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 T(z + 1) is explicitly known at the integers and half integers. The thi rd method is another recursive procedure based on the identity T(z + 1) = zT(z). The fourth method uses the Chebyshev version of a finite Fourier transform on the closed form expression for JE,T from Sec-t ion 4.1.2. This last method is essentially a numerical approximation of the integral (6.1). 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 ajt(r) = (—l)kCk(r).) (6.1) 6.1 The Lanczos Method 72 Chapter 6. Computing the Coefficients First define Tr[z) as the integral function J7, w - / > i-x2)zfEAx) dx xL appearing in equation (3.14), and notice by this same equation that Fr can also be writ ten Tr(z) = 2-V2T{z + l / 2 ) ( z + r + l / 2 ) - z - 1 / 2 e x p ( z + r + 1/2) . Now denote the 2kth Chebyshev polynomial by T2k(x) — Y^j=o CwkX2j, and from (B.2) the coefficients in the Chebyshev series are defined as 2 f1 °k{r) = ~ fE,r{X)T2k{x) n J-l dx X* The Chebyshev series coefficients are then given by "l dx fE,r(x)T2k(x) — -1 2 fl ck{r) = ~ j fE,r(x)T2k(x) x< 2 r * J-l 2 f1 = - \ JEAX) n J-l A—n J—1 dx dx 3=0 k xz 3=0 using (B.3) (6.2) W i t h a ready list of Chebyshev polynomial coefficients and precom-puted Tr values, (6.2) can be concisely expressed 2 TT Co,o C*0,2 C'2,2 " ^r(O) " ^r-(l) co(r) - c i ( r ) _ C*0,2fc C*2,2fc - ' C*2A;,2fc T r ( k ) . ( - l ) f c c f e ( r 73 Chapter 6. Computing the Coefficients In other words, Ck(r) is simply a weighted sum of J-r 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 (2n)-1/2T(z+l)(z+r+l/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 C2jt2k is required. For j = 0, C 0 ,2fc = (—l)k, while for j > 1 C ^ = ( - l ) ^ ( * ^ ) 4 « . Using this in equation (6.2) yields CQ(T) = ^2e/[n(r + 1/2)] er and The corresponding a^(r) appearing in (1.1) are then ak(r) = (—l) f cCfc(r). Al though 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 wi th 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 ce(6) using (6.2). In Table 6.1 the individual terms ar,e listed; note the order of the largest scaled sum-mands. Yet the sum of the last column is only 1.711 x 10~ 1 0 , which when multiplied by 2e6/n, yields c 6(6) = 0.00000004396. 6.2 A Horner Type Method This method is based on the simple observation that the series Sy (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) = T(z + l)(z + r + l / 2 ) - ( 2 + 1 / 2 ) e 2 + r + 1 / 2 ( 2 7 r ) - 1 / 2 , (6.3) 74 Chapter 6. Computing the Coefficients j C2j,12 0 1 .810495300735 .810495300735 1 - 7 2 .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 stated 1 1 Fr{z) = 2 a o ( r ) + *MjTi + °*<r) (,+VI 2) + from which the coefficients can be efficiently found using a type of Horner procedure1: a 0 = 2F r (0) Oi (^r(l) ap_\ 2 2 / 1 o 2 and in general, for n > 3, «n(r) = ^ ^ ( F r ( n ) - y j — - axj — - a2j — an.x (6.4) Al though 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. Fr(z) is clearly identical to Sr(z). The use of Fr(z) is meant to indicate that 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. A s such, it inherits the standard properties of this function, in particular the fundamental recursion T(z + 1) = zT(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, H k { z ) _ r(, + i)r(, + i). T(z - k + l)T(z + k + 1) z • • • (z - k + 1) (z + l)---(z + k) It is a simple exercise to deduce the following relationships: (6.5) Hk(N-l)={N k ^ 2 N + k)Hk(N), 0<k<N-l (6.6) and H^=w^vm^yw °-k-N- (67) Now oo T(z + 1) = (z + r + l/2)z+l/2e-^+r+l^V27rY^'ak(r)Hk(z) , (6.8) fc=0 and oo •zT(z) = z(z + r-l/2)z-l/2e-iz+r-l/2)V2^Y^ak{r)Hk(z-l) . k=0 B y the fundamental recursion, the right hand sides of these two equa-tions are equal, from which it follows that ez(z + r - 1/2Y-1'2 V' ak{r)Hk{z - 1) 1 = ^ : . (6.9) (z + r + l / 2 ) * + V 2 ak(r)Hk(z) For ao(r), set z = 0 in equation (6.8) to get ao(r) = y/2e/[n(r + 1/2)] t For N > 1, set z = N in (6.9). Bo th series then terminate, one at the 76 Chapter 6. Computing the Coefficients N — 1 term while the other at the iVth term. Isolating ajv(r) and simplifying using (6.6) and (6.7) then yields J V - l aN(r) = (2N)\ £ k=0 •-1/2\N~1/2 (N-k)(N + k) • + 1/2J N(N + r + 1/2) ~ ak(r) (N-k)\(N + k)' A s wi th 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 (2mkx)}kxL_00 is orthogonal not only wi th respect to in-tegration over [0,1], but also wi th respect to summation over equispaced data on this interval. Specifically, in the case of integrals ^ e - ^ { J - - „ 1 0 ) For the discrete analogue, fix an integer N > 2 and consider the values X j = j/N, j = 1 , . . . , N in [0,1]. Then for 1 < n , m < N, 1 N ( 1 - f — \ 2mnxj —2mmxj I 1 1 n m i (P, ^^\ N ^ " 1 0 else . { ° ' 3=1 K More generally, for n,m G Z , N _1 e2ninxie-2nimXj _ J 1 if W = 771 (mod TV), , , " 0 else . 3=1 V 77 Chapter 6. Computing the Coefficients Now assume that f(x) has a uniformly convergent Fourier expansion f(x) = 2\2%-oo bje2mix on [0,1], and that the values of / are known at the points Xj = j/N, j = 1,...,N. Compute approximations to the Fourier coefficients as N k=l and form the discrete Fourier series N i = i Then / = / on the equispaced data points Xj, j = 1,...,N, and / is a good approximation to / for the points in between. The larger the number TV of sampled data points becomes, the better the approx-imation. This is clear from equation (6.11) since this expression is an iV-subinterval Riemann sum approximating the integral in (6.10). In fact, it is possible to derive an explicit relationship between the coeffi-cients bj of the approximation, and the bj of the Fourier series. Since f{xn) = f{xn)i n = 1,..., TV, it follows that N oo • ' bke2nikXn = bke2nikXn , fc=l fc=—oo again for n = 1,..., N. Now multiply both sides by exp (—2mjxn), sum over n = 1,..., JV, and use (6.12) to get h = £ bk • fc=j (mod JV) For this reason, the smoother the function / being approximated, the more rapidly its Fourier coefficients bj decrease to zero, and hence the better the one term approximation bj ~ 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 TV zeros x\,..., at xk = cos I — Define the coefficients Cn by N TV 2 Cn = — ^ / ( x f c ) T n ( x f c ) fc=l 4 t / ( c o < ^ ) ) c o s ( ^ ) ) fc=l Then the approximation N f{x)^Y'^T^x) n=0 is exact on the TV zeros of TN(X). Larger TV values produce more accu-rate estimates of cn and hence f(x). This method can be applied to estimate the coefficients in (1.1.) using the explicit form of fE,r(%) from equation (4.4): 2 N an = (-l)nCn = ( - 1 ) " - Y fEAxk)T2n(xk) . fe=l This calculation can be economized since fE,r{xk) a n d Tin{xk) are even wi th respect to xk, hence, taking TV even, ' - ' r i v * W ( / 7 r ( f c - l / 2 ) \ \ / 2 m r ( f c - l / 2 ) A" = (-1) N H" [COS { TV J J C ° S { TV Furthermore, the TV values / ^ ( x i ) , . . . , fE,r{xN) need only be com-puted once since they do not depend on n, the index of the coefficient being computed. Selecting TV = 20 was sufficient to reproduce do,... ,0,5 to six dec-imals as given in Lanczos' original paper for r = 1,1.5,2,3. In fact, TV — 20 is really only needed for r = 1, as the error \cn[r) — cn\ appears 79 Chapter 6. Computing the Coefficients to drop drastically as r increases. For example, wi th r = 3 and N = 20, 10 decimal precision is achieved for the first 6 coefficients. Wha t 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-\(y) and Wo(y) near their branch points y = 0 and y = —1/e. In [5] the authors discuss the subtleties involved wi th evaluation near these points. The numerical experiments noted here were carried out using Maple 8 which has bui l t - 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 a 0 / 2 +1.4598430249 +2.0844142416 +3.0738046712 +7.0616588080 a i -0.4606423129 -1.0846349295 -2.1123757377 -6.5993579389 « 2 +0.0010544242 +0.0001206982 +0.0386211602 +0.5396522297 -0.0003384921 +0.0001145664 -0.0000510050 -0.0019519669 04 +0.0001175425 -0.0000176145 +0.0000004776 -0.0000013258 a 5 -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 r N T{z + 1) w (z + r + l / 2 ) 2 + 1 / V ( 2 + r + 1 / 2 ) ^ r " (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 Sr(z) = 1 . . , \ z , \ z(z — 1) 2 " ° ( r ) + a ' ( r )7TT + a 2 ( r ) ( , + i ) ( , + 2 ) + = ^2 ak(r)Hk(z) . k=0 If this is terminated after the TV'th term, the goal is to express the resulting series Sr^(z) — ao(r)/2 + X^j t l i ak{r)Hk{z) in the form SrAz) = b0(r) + Y ^ l • Observe that for k > 1, bk(r) is the residue of Sr>N(z) at z = —k, while, upon taking z —> oo, JV bo{r) = a0(r)/2 + Yak{r) • k=i -Now Res [ t f f c ( z ) ] 2 = , = T ; (fc-j)'-[0-UHa ' ^ 10 else , so that To determine 6j(r) of equation (6.13), simply sum up the coefficients of l/(z + j) in SrtN(z): K—j 81 Chapter 6. Computing the Coefficients 6.7 Matrix Representation A s remarked in Section 6.1, the calculation of the coefficients can be concisely expressed in matrix form, which, along wi th the partial frac-t ion 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 ao ( r ) , . . . ,apf(r), and let b(r) be the column vector of partial fraction coefficients bo(r),..., 'jjv(r). Wri t ing C as the matrix of Chebyshev coefficients C ^ j , 0 < i, j < N, C has entries Cij (i+j)(2j)!(i-j)! Define the partial fraction decomposition matrix to be B , which has entries ' 0 if i > j, 1/2 iii = j = 0, 1 if i = 0, j > 0, as given in Section 6.6. Denote by F r the column vector of scaled values \/2Jrr(j)e~r~1/2, 0 < j •< N, where as in Section 6.1, Tr(z) = 2~1/2F{z + l/2)(z + r + 1 /2 ) - 2 " 1 / 2 exp (z + r + 1/2) . W i t h this notation, the Fourier coefficients appearing in (1.1) can be concisely expressed a(r) = ^ C F r 7T while the coefficients for the partial fraction decomposition are simply x / 2 e r + 1 / 2 b(r) = — B C F r . 7T 82 Chapter 6. Computing the Coefficients , The only floating point calculations are in the term ( y 2 e r + 1 / 2 / 7 r ) F r , so that if multiple sets of coefficients are required for different values of r , only this term need be recomputed. St i l l further efficiency can be realized by observing that the leading e r + 1 / 2 term of the b(r ) coefficients cancels once inserted into (6.13). Let t ing u(z) denote the vector of rational functions 1, ( z + 1 ) - 1 , . . . , (z+ l)~k, equation (6.13) becomes simply e f z + r + 1/2 ) 2+1/2 r(z + i) ^ 2 u T ( z ) B C F r . (6.14) 83 Chapter 7 Lanczos' Limit Formula The final result of Lanczos' paper [14] is the beautiful l imit formula, stated here as T h e o r e m 7.1 . For Re(z) > 0, oo T(z + 1) = 2 l im rz S"'(-l)ke-k2/rHk(z) . (7.1) 1 KX ' ' k=0 This result 1 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 vr. In this chapter the validity of (7.1) is proved in detail for Re(^) > 0. Note that the formula is misstated in the review [29]. There the term rz is incorrectly stated as r 2 . 84 Chapter 7. Lanczos' Limit Formula 7.1 Some Motivation A s 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 d0 cos 2 2 0 •TT /2 logv _ ^ + r + i / 2 y + 1 / 2 y 7 r / 'z + r + l / 2 \ ^ r / 2 c o s ^ 2 ( v / e r s i n ^ ^ ' - T T / 2 logv A rescaling of r to er produces the r 2 term of equation (7.1) and elim-inates the exponential term in the limit, thus T ( z + 1 / 2 ) = 2r« + « + f'2 c o s , , r^v/er^e d e \ er ) J_^/2 logv = 2r*(Z + 6 r e r 1 / 2 y f cos2z0r1/2gr(0)d0szy, ^ (7-2) where Now gr(0) converges point-wise to.zero rapidly on [—7r/2,7r/2) wi th increasing r , but equals one at 0 = TT/2. A s a result, as r increases, (r/ny/2gr(0) takes on the appearance of a the Gaussian distribution ( r / T r ) 1 / 2 exp [-r(0 - TT/2)2]. Plots of exp [-r(0 - TT/2)2] - gr{0) for in-creasing values of r show a convincing fit; see Figure 7.1. Furthermore, 4= / cos ^ y / V ^ / 2 ' 2 dt = (-l)ke-k2/r . V71" J-oo so that if gr(0) is replaced by exp [—r(8 — 7r/2)2] and the bounds of integration [—7r/2,7r/2] are extended to ( — 0 0 , 0 0 ) , the integral gives rise to the terms e~k lr upon writing cos 2 2 0 as its Fourier series. 85 Chapter 7. Lanczos' Limit Formula These ideas lead to the following outline of the proof of equa-t ion (7.1). For large r and fixed z, T(z + 1/2) = 2rz (^  + e r + 1 / 2V+ 1 / 2 cos2* flr1'2 M TT/2 2rz C ' \ o s 2 z Q r l l 2 e - r ^ W 2 d6 -TT/2 rz I \ o s 2 z 6 r l l 2 e - r ^ W 2 d6 I. 86 Chapter 7. Lanczos' Limit Formula \ 0 / 2 r (z + i/2) ' V^F r (z + i) # f c(z) cos (2fc0) r1 / 2 e-'WV d6 .,r(* + l /2) r (z + 1 ) 0 0 = 2r X)'(-l) f cc- f c a / r/r f c(z) , 0 where the « symbol becomes equality as r —> 00 . 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 wi l l be justified by a series of smaller lemmas. Begin wi th a definition: D e f i n i t i o n 7 .1 . For r > 0, define Clearly 5r —> 0 as r —> 00 . It wi l l also be convenient to introduce some notation for the various intervals of integration: D e f i n i t i o n 7.2. For Sr < 7r, define the intervals 1/4 (7.3) 7T 7T and 2 ' 2 The first lemma is straightforward: 87 Chapter 7. Lanczos' Limit Formula L e m m a 7.1. The function 0 if 9 = 0, h(9) = , • „ v > ^ ml eise logv is continuous on [—TT/2, TT/2], and is thus bounded. P r o o f of L e m m a 7.1: O n (-TT/2, TT/2], h(6) = fer(9)v-er2-1/2, which is continuous by Theorem 4.3. A t 9 = —TT/2, l im h(9) = 0 = hi-TT/2) . The next lemma shows that the integral in equation (7.2) is con-centrated on the small interval Ir containing 7r/2: L e m m a 7.2. Suppose Re(z) > 0, and let Jjr ^eJ logv Then l im rz T\ (r, z) = 0 . 1—>oo P r o o f of L e m m a 7.2: First , recall the calculus result s in 2 x l im —r— = 1 , so that for \x\ sufficiently small, s i n 2 x > \x2 . (7.4) 88 Chapter 7. Lanczos' Limit Formula Now bound X\ (r, z) as follows: 'v\er sind 2 / c o s 2 2 ^ 1 / 2 Qf logv de • j f cosher1'2 Q ) ' de for some C\ > 0 by Lemma 7.1 ^ ^ a + l / 2 (Vi\Tv/2-Sry\ Cr . . . . < G i r ^ ' I I | j r | since v is increasing on J r < C 2 r C T + V 2 e e r [ l o g ^ ( 7 r / 2 - 5 r ) - l ] = c 2 r < T + 1 / 2 e - e r c o s 2 ( 7 r / 2 - ' 5 r ) / , ; by the implicit definition of v{6) < C 2 r C T + 1 / 2 g - e r c o s 2 ( 7 r / 2 - 5 r ) / e = C 7 ' C r + : ! - / 2 g _ ' " s i n 2 ^r-< C 2 r a + l ' 2 e - r S - / 2 for 5 r sufficiently small, by (7.4) = C2ra+1/2e~{1/2)^r/logr 0 as r —>• oo The next lemma justifies the replacement of gr(0) by exp [—r(8 — 7r on J r Lemma 7.3. Suppose Re(z) > 0, and let / \ f or „ i l o I/"v\er sm# l2(r,z) = / cos eT (-) JIr \_\eJ logv _ e-r{B-Kllf de 89 Chapter 7. Lanczos' Limit Formula Then l im rzZiir, z) = 0 P r o o f o f L e m m a 7.3: Begin wi th fv\er s in0 rz I cos 2 2 Or1'2 IT ^vy a m p _ ^ _ r ( 9 _ 7 r / 2 ) 2 - e / log v dO <ra f Jlr cos2" O r 1 ' 2 e - ^ ' 2 ^ V e / log v de Lett ing Xr denote the indicator function on Ir, it is enough to show that as r —> oo (i) uniformly, and (ii) that v\er s in0 0 e ) log v er ( 0 - V 2 ) 2 _ l Xr —>0 Si, r° / cos2° 6 r ^ e ' ^ - ^ d O remains bounded. For (i), since l ime_ + 7 r /2- s i n 0 / l o g v = 1, it is enough to show that (§)' er(e-7r/2)2 _ jj = r e r( lo g U - l )+r(e -7r /2) 2 _ 1 ' uniformly, which boils down to showing that [er(logv - 1) + r(0 - TT/2) 2] ^ — > 0 uniformly. Expanding e(logv — 1) + (0 — 7r /2) 2 about 0 = 7r/2 gives e ( l o g v - l ) + ( 0 - 7 r / 2 ) 2 = (l~l)V- *W + ^ ( " 4 - f + f ) (« - + O ((« - ,/ 90 Chapter 7. Lanczos' Limit Formula Now the coefficient (1/3 — 1/e) < 0, so that for |0 — TT/2\ sufficiently small, - ( 0 - T T /2 ) 4 < e{logv - 1) + (0 - T T /2) 2 < 0 , so that, for r sufficiently large, - ( 0 - T T /2 ) 4 < [e(logv - 1) + (0 - TT/2) 2] Xr < 0 , whence - r ( 0 - T T /2 ) 4 < [er(logv - 1) + r(0 - TT/2) 2] Xr < 0 . The term (0 — 7r/2) is at most 5 r in magnitude on Ir, so - r 5 4 < [er(logv - 1) + r(0 - TT/2) 2] Xr < 0 , and hence "r (rloW) " M ^ - ^ + ^ - v r / S ) 2 ] Xr < 0 . Now let r —> oo to get [er(log u - 1) + r(0 - TT/2) 2] Xr — • 0 uniforly. This proves (i). For (ii), first observe that about 0 = 7r/2, C O S 2 c t 0 < ( 0 - T T / 2 ) 2 < t . Thus rC T / c o s 2 C T 0 r 1 / 2 e - ^ 9 - ^ 2 d 0 J l r <r° f ( e - T T ^ r ^ e - ^ - ^ d d J l r poo <r° / ( 0 - 7 r / 2 ) 2 C T r 1 / 2 e - ^ - 7 r / 2 ) 2 d 0 . 91 Chapter 7. Lanczos' Limit Formula Now make the change of variables u = r(6 — iv/2)2, du = 2r(9 — n/2)dd. This gives /•CO • rC T / ( d - n ^ r ^ e - ^ - ^ d e . Jn/2 = 7i / W-X'2e-Udu 2 Jo That is, r° [ c o s 2 c r 6 ' r 1 / 2 e - r ( e - 7 r / 2 ) 2 dO < \v{o + 1/2) . Jlr 2 This proves (ii) and completes the proof. L e m m a 7.4. Suppose Re(z) > 0, and let J 3 (r, z) = /" cos 2 2 0 r ^ e ^ " 7 ^ 2 d8 . Then l im rz Tz (r, z) = 0 . 92 Chapter 7. Lanczos' Limit Formula Proof of Lemma 7.4: rz f COB2'Br1'2 e-^'^dB <r° [ rl'2e-^-^2de poo = ra e-" 2 du , J y/r6r this last line a result of the change of variables u = y/r(9 — TT/2). This can be expressed in terms of the error function erf (x) = —= / e -* 2 dt V 7 1 " Jo as x —> oo . Picking up where we left off then we have r° e~u2 du = r a ^ - [l - erf(Vr"5r)] 2 ^/r8r = I r ^ - l / ^ l o g r ) l / 4 e - ( r / l o g r ) V 2 —> 0 as r —• oo , completing the proof. 93 Chapter 7. Lanczos' Limit Formula L e m m a 7.5. For integer k, 1 r°° - 7 = / cos(2key^e-r^-7C^2d9= (-l)fce-fc2/r . (7.5) V7T 7_oo P r o o f o f L e m m a 7.5: The integral in (7.5) is the real part of i r°° J_ / r l/ 2 e i 2 f c f l - r ( f l - V 2 ) » dy V 7 1 " J-oo = (_i)ke-k2/r_}_ r r i / 2 e - r ( e - ^ t - ) 2 de u p o n completing the square 7-00 - ( _ i ) ^ e - f c 2 / ' - _ L / r 1 / 2 e~ r 6 d6 wi th a simple contour shift y/n 7-oo = ( - l ) f c e - f e 2 / r . We are now in a position to prove (7.1). P r o o f o f T h e o r e m 7.1: Recall equation (7.2) which is true for al l r > 0: r ( z + 1 / 2 ) = 2r« (* + ~ + l / 2 \ 1^ c o s- 1 r^lmi^l i0 . V e r J J-n/2 ^gV From the definitions of T\(r, z), J 2 ( r , z) and Tz(r, z) in Lemmas 7.2, 7.3 and 7.4, respectively, v ( - + e r + l/2^ ^ ^ ^ r v , ( / e r d a > ^ T/2 logV = 2r ^z + er + l/2y+1'* J*'' z + er + l/2\z+1/2 er x 1 z*00 I i ( r , ^) + J 2 ( r , 2) - 2J 3 ( r , z) + ^ / cos 2 2 flr1/2 e ^ * - ^ ) 2 d 6 , ^ J - 0 0 94 Chapter 7. Lanczos' Limit Formula Now let r —> oo and apply Lemmas 7.2, 7.3 and 7.4 to get 1 f°° Y{z + 1/2) = 2 l im r z - \ cos2z Or1'2 e " ^ - ^ ) 2 d 9 _ »--°o 2 J _ 0 0 Finally, replace cos 2 2 0 wi th its Fourier series (see Section 3.3) OO c o s 2 2 0 = r t + y ^ ]Hk{z) cos {2k6) and integrate term by term using Lemma 7.5 to get T(z + 1/2) = 2 l im r » r ( * + * f f f ^ - l ) f c e - f c 2 / r # f c ( * ) . r^oo 1 (z + 1) ^ Canceling the T(z + 1/2) terms and moving T(z + 1) to the left hand side completes the proof. 7.3 Additional Remarks In his original paper [14], Lanczos claims that the l imit in Theorem 7.1 converges for al l 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 l imit formula on the l imit ing behaviour of the individ-ual coefficients ak(r), not on the l imit ing behaviour of the integral (7.2). Specifically, the comparison appears to be 0 0 T(z + 1) = v 7 ^ (z + r + l / 2 ) 2 + 1 / 2 e - ( z + r + V 2 ) ^ ak(r)Hk{z) fc=0 V ^ fc=0 v CO « 2rzY^(-^)ke~kVrHk(z) . o 95 Chapter 7. Lanczos' Limit Formula This suggests that for large r, fnrak(er) 2 e e or, by rescaling r , Q(k,r) = ( _ l ) V " e f c > J - M r ) ) - 1 ~ 1 . Tfr Plots of logQ(/c, 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 J—•—T~ ~1~rT-i L fc 24 Figure 7.2: logQ(/c,10), k = 0 , . . . , 24 96 Chapter 7. Lanczos' Limit Formula -35' • 0 k 24 Figure 7.3: log Q(k, 20), k = 0 , . . . , 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 T(z + 1), the series (1.2) is truncated and an estimate of the ta i l of the series eTtn(z) = YlT=n+i ak(r)Hk(z) is therefore required. The error function er,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-t ion of some new notation and an important theorem, Lanczos' own error estimates and those of Luke [17] are examined; The final sec-t ion deals wi th numerical investigations motivated by the observation that the error function er,n(z) is the tai l 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 wi l l be shown to occur on the imaginary axis. Furthermore, the value of this maximum can be easily estimated empirically wi th the mapping t —» it/'(1 — t) of the interval [0,1) onto the positive imaginary axis. Finally, the l imit ing behaviour of er,n(z) as z —> oo 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 | e 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 er>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. A t 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 wi l l be in the context of the gamma function. A s such, let G(l + z) denote an approximation to F ( l + z), and let e a denote the difference between the approximation and the actual value of the function: This is often termed the "absolute error", even without taking the modulus of ea. The relative error is here denoted ep(z), and it is defined according to the commonly accepted notion of (actual-estimate)/actual: From the definition of ep(z), and assuming \ep(z)\ < 1, r(l-f-z)'may be expressed ea(z) = T(l + z)-G{l + z) . T{1 + z) - G(l + z) T(l + z) (8.1) T(l + z) • r(l + z) G(l + z) 1 - ep(z) G(l + z)[l + ep(z) + ep(z)2 + ---) (8.3) G(l + z)(l + e„(z)) . (8.4) 99 Chapter 8. Implementation & Error Discussion It is common in the literature to encounter expressions of the type T ( l + z) = G(l + z)(l + e(z)), where e(z) is termed the relative error. This is not, strictly speaking, relative error in the sense of (8.1), but comparison wi th 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 wi thin 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 X\ = 123 and x2 = 0.123 are computed, both wi th a relative error of ep = 0.1. Then the actual errors are x / c e p / ( l — e p), k = 1,2, which are approximately 14 and 0.014, respectively. This means that for both x\ and x2, the first non-zero digit is good while the second may be in question. B y contrast, if X\ 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 T(z + 1) = (z + r + i / 2 ) 2 + V 2 e - ( 2 + r + i / 2 ) ^ [ S r t n { z ) + € ^ z ) ] ( g 5 ) where Sr,n(z) is the series (1.2) truncated after a finite number of terms, and he calls trjn(z) the relative error. It turns out that Sr,n(z) f» 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 Spouge makes this same observation in [27] and shows that for Re(z) > 0, the denominator is bounded below by T(z + l)(z + r + l / 2 ) - ( ^ + 1 / 2 ) ^ + r + i / 2 ( 2 7 r ) - i / 2 • \T(z + l)(z + r + l / 2 ) - ( 2 + 1 / 2 ) e 2 + r + 1 / 2 ( 2 7 r ) - 1 / 2 | > 100 Chapter 8. Implementation <fe Error Discussion where upon selecting r = z = 0 the bound in this inequality is sharp. Thus /7T\ V 2 \eP\ < [-) \er,n(z)\ , . (8.6) and since \J(ir/e) = 1.075, bounding relative error in the Lanczos sense gives essentially the same bound for ep, the standard notion of relative error. For the purposes of this study z wi l l be restricted to Re(z) > 0 and as such the Lanczos notion of relative error wi l l be used. It should be pointed out, however, that wi th er,n(z) something is lost if z is moved to the left of the imaginary 1 axis. A s we saw in Chapter 5, the formula (1.1) extends to the region Re(z) > —r not including the negative integers. Lanczos' relative error er,n(z) has poles at the negative integers in this case, and becomes unbounded. ep(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: r(i + *)r(i-*) = ^ - . sin TTZ Already this says that methods used to estimate T(l + z) need only be valid on Re(z) > 0. Now suppose r(l + z) is computed with a relative error of ep, that is, G(l + z) T(l + z) 1 - € p Then by the reflection formula, TTZ 1 — en T(l-z) = sin7rz G(l + z) TTZ 1 siXiTTzG{\ + z) ( l + e p + e 2 + . . . ) TTZ 1 1 sin7T2; G(l + z) 1 + ep 101 Chapter 8. Implementation & Error Discussion In other words, to compute T ( l — z) to a relative error of —ep, it is sufficient to compute r(l + z) to a relative error of ep, and vice versa. Thus when using a particular method to compute F(l+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 wi l l 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): HQ(z) = 1 , and for k > 1, . T(z + l)T(z + l) Hk(z) T(z-k + l)T(z + k + l) z---(z-k + l) (z + l)---(z + k) Observe that Hk(z) is a meromorphic function wi th simple poles at z = —1,..., — k and simple zeros at z = 0,..., k — 1. Also, Hk{z) 1 as \z\ —> oo. 2. Write (1.1) as T(z + 1) = V27 (z + r + l/2)z+l'2 e - ( ^ + i / 2 ) [ S r ^ z ) + where n Sr,n(Z) = X] ak(r)Hk(z) k=0 and oo er,n(-z) = Yl ak(r)Hk(z) . 102 Chapter 8. Implementation & Error Discussion 3. Let Fr(z) = T(z + l)(z + r + l/2)^z+1^ez+r+1/2(2n)-1/2 . B y Stirling's formula, Fr(z) —• 1 as \z\ —> oo in any half plane Re(z) > o. Note that £r,n(z) = Fr{z) - Sr,n(z), (8.7) and that again in any half plane Re(z) > cr, in particular for z in the right half plane Re(z) > 0. In view of this last l imit , the error at infinity is denoted n 4. Let n rir,n{6) = fEM - ^ ' a f e ( r ) c o s ( 2 ^ ) , fc=0 the error in the (n + l )- term Fourier series approximation of fE,r(6), where JE,r{Q) n a s t n e meaning of (3.7). Then r7 r > n(—7r/2) is the Fourier error at 0 = —7r/2, which is, incidentally, equal to e r , n ( - l / 2 ) . From equation (8.7) we state L e m m a 8.1. Forr > 0, ern(z) is analytic as a function ofz in Re(z) > - 1 / 2 . P r o o f o f L e m m a 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 e r > n (z) . • The following lemma is a consequence of the maximum modulus principle: 103 Chapter 8. Implementation &, Error Discussion L e m m a 8.2. Suppose f is analytic in the right half plane f2 = {Re(z) > o} and has the property l im |/(,z)| = C < oo . (8.8) \z\—>oo Let B = supz€Q \f(z)\. Then B = s u p R e ( 2 ) = C T \f(z)\. P r o o f o f L e m m a 8.2: Since / is analytic on fl, 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 = l im^oo \f(o + it)\; otherwise (ii) B > C. For this second case, for any compact K C O, denote by MK the maximum of | / | on K, and consider the sequence of closed right-half semi-disks {fij}£Li centered at the origin, each £lj having radius j and diameter on the line Re(z) = o. B y the maximum modulus principle, Mn. = Man,- , and MQQJ increases monotonically to B. Write the boundary of Qj as the union of two closed segments dttj = Ij UAj, where Ij is the diameter and Aj the semicircular arc of dflj. Now l im Man, = B , 3-*oo so that l im m a x { M / . , M A . } = B , but l im MA.=C <B . Thus l im Mi. = B J-+0O completing the proof. 104 Chapter 8. Implementation & Error Discussion T h e o r e m 8.1. For r > 0 and Re(z) > 0, the supremum of \er^n(z)\ in this region occurs on the line z = it, possibly at z = zoo. P r o o f o f T h e o r e m 8.1: Recall that er>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 k r , n ( ^ ) | 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 friO)1 at 9 = —TT/2, that is by forming the alternate sum 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 rj. Then a simple estimation shows that the influence of this error on the integral transform (16) (for values of z which stay wi thin the right complex half plane), cannot be greater than (n/2)r]. Thus we can give a definite error bound for the approximation obtained. 1 In his paper, Lanczos uses 7 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. n 105 Chapter 8. Implementation &. Error Discussion The "integral transform (16)" referred to here is the integral in (3.6): r+n/2 / cos2z (6) fr (9) dO , and Lanczos' comments relate the error in truncating the series in (1.1) after the n t h term, e r , n (z), to the error T]r^n(9) in the (n+l ) - t e rm Fourier series approximation of fE,r{9)- Based apparently on this observation, he gives uniform error bounds for Re(z) > 0 as listed in Table 8.1. n r 1 1 0.001 1 1.5 0.00024 2 2 5.1 x 10~ 5 3 2 1.5 x 10" 6 3 3 1.4 x 10" s 4 4 5 x 10~ 8 6 5 2 x 10~ 1 U Table 8.1: Lanczos' uniform error bounds The steps connecting Lanczos' observation about r]r,n(9) 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(2) | < 5 x 10~ 8 . It is a simple matter to compute directly l im |c4 > 4 ( t t ) | = l - ^ T a t ( 4 ) t—>oo 2 —' fc=l = 5.3 x 1 0 - 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{9) had a finite jump discontinuity at 9 = —TT/2. In this case, however, fE,r(9) is continuous at the endpoints 106 Chapter 8. Implementation &; Error Discussion n r 9r,n 1 1 -it/2 1 1.5 0 2 2 -TT/2 3 2 « -0 .72 3 3 -TT/2 4 4 0 6 5 - T T / 2 Table 8.2: Location of maximum of r)r^n(9) n r 1 ,oo 1 1 8.0 x 10~ 4 1 1.5 2.2 x 10" 4 2 2 5.0 x 10" 5 3 2 9.1 x 10" 7 3 3 1.1 x 10~ 6 4 4 5.3 x 1 0 - 8 6 5 1.9 x 1 0 " 1 0 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 6rtn denote the least value of 9 in [—7r/2,7r/2] where the maximum of r]^n(9) occurs, we find (empirically) the values listed in Ta-ble 8.2. Indeed, it is not difficult to construct examples in which rjr>n(9) is zero at 9 = —TT/2. Take for example n = 1 and r = 1.500773 In this case the maximum of \rjr,n(9)\ is approximately 0.00024 and it occurs at 9 = 0. One may surmise that Lanczos' bounds are, for each choice of r and n, simply the error at infinity, e^°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 |e r i n(2)| does occur as \z\ —> oo. This is not the case in general though, as wi l l 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. Empir ica l 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 7/V,n(#) near singularities of /e, r(6'), and he uses this notion to bound e r , n(z) for o = Re(z) > 0. Luke's bound is reproduced here. From the derivation in Section 3, oo ennO) = ak(r)Hk(z) k=n+l i r(z + i) i r(z + i) v^r(2 + i/2) Yl °fc(r) / t=n+l J ~ « / 2 /•TT /2 / cos 2 2 dd J-ir/2 cos2z6cos2k0d0 Thus, assuming |// r,n(^)| is maximized at 0 = — TT/2, \^r,n(z)\ — i— 'TT r(z + i) r(z + i/2) / C O S 2 < T % r , n ( - 7 T / 2 ) M 0 J-TT/2 1 7^ r)r,n(-n/2)r(z + l) T(z +1/2) r(q+-l)r(q + l/2) T{o + l)T(o + l) ^n(-jr/2)T(z + l) T(z +1/2) r(q + l/2) r(a + l) Unfortunately, this argument breaks down when one considers that Vr,n(—'K/2) has zeros as a function of r , as was pointed out i n 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 T(z+1) on the right half plane. Furthermore, this error bound is not uniform since it grows without bound as \z\ —• oo along lines parallel to the imaginary axis. This subtle role of r on the error bounds motivates a closer exami-nation of ertTl(z) as a function of r. 8.5 Numerical Investigations &; the Zeros of e°° 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 er<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 MQ of |e r , n (z) | in the right half plane occurs on the imaginary axis, possibly at z = ±ioc. 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, e ^ 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 k. Error Discussion case ( i ) : Mn = 0 This case is impossible, for otherwise the result would be a closed form expression T(z + l) = V2^(z + n + l/2)z+l'2 e^z+r'+1^Srt,n(z) which could be continued analytically to the slit plane C \ (—oo, — r* — 1/2), but which would remain bounded near — (n + 1), — (n + 2 ) , . . . , contrary to the unboundedness of T at the negative integers. case ( i i ) : M n > 0 In this case the maximum of | e r „ ) n ( 2 ) | occurs on a bounded segment of the imaginary axis; the question remains: how to locate it? Exper i -mentally, for each n = 0, . . . . , 10, zeros r* were found and |er,)Tl(2;)| 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 kr„,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/(l — t), 0 < t < 1, which maps the unit interval onto the positive imaginary axis, the rational functions Hk{z) take the form Hk(z(t)) fc-i =n 3=0 z(t)+j + l fc-1 j i "/a -i t / ( l - t ) - j t/(l-t)+J + l TT t(i + j) - j 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 er>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 er.,n(>(*)) ~ Y ak(r*)Hk(z(t)) k=n+l for M not too large. The problem of bounding | e r „ ) n ( z ) | in the right half plane reduces to that of estimating the maximum of 2~2kLn+i ak{r*)Hk(z(t)) for t e [0,1). W i t h this strategy, the existence of zeros of e™n is no longer an issue; for any r > 0, simply examine the finite series approximation of \er}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 produced the smallest uniform error bound for \er^n(z(t))\. 8 .5 .2 Z e r o s o f c°° The question remains, for n fixed, does e~n have zeros as a function of r? The answer was found to be yes in all cases examined. Refer to F ig -ure 8.1 for plots of [ e^J 1 / 5 , n = 0,1, 2, 3, showing the location of zeros for the first few cases.2 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 efn n = 0 , . . . , 60. The results of this investigation are listed in Tables C . l and C.2 of Appendix C . From these tables, e^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 I m p r o v e d L a n c z o s F o r m u l a A s 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]. The functions 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. Ill Chapter 8. Implementation & Error Discussion That is, r(z+l) = V^(2+r+ l /2 ) 2 + 1 / 2 e - ( 2 + r + 1 / 2 ) ]P ak(r)Hk(z) + e r , 6 ( » k=0 W i t h the choice of r = 5, Lanczos gives a uniform error bound of |e r , 6 (z) | < 2 x 1 0 " 1 0 . Twelve zeros of e£°6 were found; call these r o , r i , . . . , r n . For each of these zeros, using the estimation techniques of Section 8.5.1, the maximum Mrjt& of | e r j i 6 ( z t / ( l — t))\ was estimated along wi th its loca-t ion on the imaginary axis using first fifteen terms of | e r j i 6 ( i t / ( l — 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 x 10~ 1 2 , nearly 1/100 of Lanczos' estimate. j t/(l-t) 0 -0.117620 0.565599 4.71 x 10" 4 1 0.684391 1.525502 2.75 x 10" 6 2 1.450013 2.388058 8.88 x 10" 8 3 2.182290 3.172714 6.78 x 1 0 - 9 4 2.883225 3.887447 9.30 x 1 0 " 1 0 5 3.553321 4.538907 1.99 x 1 0 - 1 0 6 4.191832 5.133674 6.07 x 1 0 ~ n 7 4.796781 5.678871 2.49 x 1 0 - 1 1 8 5.364813 6.183352 1.30 x 1 0 ~ n 9 5.891184 6.660957 8.02 x 1 0 " 1 2 10 6.372580 7.137483 5.29 x 1 0 " 1 2 11 6.779506 7.883760 2.72 x 1 0 " 1 2 Table 8.4: Approximate Zeros of e£ Plots of log 16^.^(^/(1 — t))\ for j = 0 , 1 , . . . , 11 are shown in F ig -ure 8.2. Observe that the maximum of the curve associated wi th the largest zero r n gives the smallest maximum error of approximately exp( -26 .6 ) . To get an idea of how the maximum error along the imaginary axis varies wi th r , refer to Figure 8.3 in which log |e r ) 4 ( i£ / (1 — t))\ is plotted f o r 0 < t < l , 0 < r < 5 . The ends of the deep troughs near t = 1 correspond to the zeros of e ^ 3 . 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 For graphing 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 8.5.4 Optimal Formulas In practice, for each value of n examined, the largest zero of was found to produce the smallest maximum error (among all zeros), and bounds much improved on Lanczos' estimates. Tables C . l and C.2 of Appendix C lists, for n = 0 , . . . , 60, the largest zero of e£°n, denoted r (n) , the corresponding estimates M r ( n ) ) n ) 5 and M r ( n ) n ) i 5 of the maxi-mum of | e r ( n ) n (z t / (1 — t)) | using five and then fifteen terms of the sum, respectively, and finally, the location of the maximum of |e r ( n ) ) T l (zi) | on the imaginary axis . 4 . From this data it appears that five terms of the series |e r(n),n(^/(l — ^))| suffice to estimate the maximum error in al l 4 O f particular 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 k. Error Discussion - 2 2 Figure 8.3: log | e M ( i t / ( l - *))|, 0 < t < 1, 0 < r < 5 but the first few cases. Observe that the n = 10, r = 10.900511 row of Table C . l is sufficient to guarantee 16 digit floating point accuracy in the right-half plane. In this case, resolving the series 5^(10) ,10 into partial fractions and rescaling (7 = 1-5): z\ = (z + 2 ) 2 + 1 / 2 e - ( 2 + 2 ) v / 2 ^ (o.999779 + -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 7 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 e^. 115 Chapter 8. Implementation & Error Discussion the coefficients as in equation (6.14) yields the approximation 10 r ( z + 1 ) K 2 ^ ( i ± l M ± l ^ 2 + 1 / 2 d0 + J2 fc=i z + k where the coefficients dk are as given in Table 8.5. k dk 0 +2.48574089138753565546 x 10~ 5 1 +1.05142378581721974210 x 10° 2 -3.45687097222016235469 x 10° 3 +4.51227709466894823700 x 10 u 4 -2.98285225323576655721 x 10 u 5 +1.05639711577126713077 x 10 u 6 -1.95428773191645869583 x 10" 1 7 +1.70970543404441224307 x 10~ 2 8 -5.71926117404305781283 x 10" 4 9 +4.63399473359905636708 x 10~ 6 10 -2.71994908488607703910 x 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 e£°n have? What is the largest zero? (ii) For fixed n, why does the largest zero of seem to give rise to the smallest uniform bound on |e r i r a(z)|? (iii) For fixed n and r, is it possible to analytically determine the maximum of \er<n(it/(l — t))\ for t G [0,1)? (iv) Perhaps most importantly, given e > 0, what is the best choice of r and n to estimate T wi th a relative error of at most e? 116 Chapter 8. Implementation & Error Discussion It was observed early in the investigation that, not only does e£°n have zeros as a function of r , but so too do the individual ak{r). This is not surprising, since oo fc=n+l assuming the coefficients are rapidly decreasing. However, the relation-ship between e^n and an+i(r) is eye opening; refer to Figures 8.4 and 8.5 for plots of [e^g]1/7 and [a7{r)]ll7, respectively. Bo th functions are then plotted in Figure 8.6 and the plots are virtually indistinguishable. ,i : : h i; li 1' i ! i i 1H j M M \i i W \ \ M M \ \ M \\ i i f J • i ' N! 1 II i i >r~~ ""!! i • ,** ;> ' 5 p f : : i i "J 1 I jl d ' 1 i iM l Hi ! t-j i; 1 j . \}< | i ...|! '-I ; ! • • : ; I 1 TTTT TiTj i ' ! •! . . L 'l 1 J i £ | • .L ' „ •; Figure 8.4: [e^] 1 / 7 , - 0 . 2 < r < 7 Numerical checks show that the zeros of e£°n and a n + i ( r ) nearly co-incide in the cases examined. Furthermore, at r(n) equal the largest zero of efn, an+i(r(n)) appears to be a very good estimate for the max-imum of \ertn(it/(l — t))\ on [0,1). This observation seems in-keeping wi th the rule of thumb for asymptotic series wi th 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 j . 4 • i \ i ' "' ' • : 1 : r ; l f | J ^ !• !' 1 K' • •• >a •* I : : : [a ! . f 1 = i j(rjji/7 i l l ! . !!! ... 1 t j l J-, S III!TTt T T T ij .. i1 1 ' |5 "Hi f1 i ! i ,j = i i j ) j i • • • j 1 l : ! L !• ! •1 1 ' )' i ' . |« i h ' "•{ } . . . . . Il ""I • Ill i ! : 1 < ; jl : !; ^ •; . ! ' " 1 r, ; i • |IU ' Ijtttfitti I 1 J Figure 8.5: M r ) ] 1 / 7 , - 0 . 2 < r < 7 li \ J i i l M , 11 ; i (i A t : : | : : : i '1' 1 rY |i i||„ I; • \ i • • it 1 1 s s : : 7 : ( " • ' |l)ij.)i;..|4| i i lp i ' 1 i i li5 1 'll P j S i l l JL . i l t]"rr i i \ -I r M I j : i* , * ' i Mi n i i. L r , i' j .1... ! i '•' 1. .1 i i ; i !: li' Figure 8.6: [e™6}1/7 and K ( r ) ] 1 / 7 , - 0 . 2 < r < 7 listed are r(n) , the estimated uniform error bound M r ( n ) ) 7 l ) i 5 , and first two terms omitted from the series, o n + i ( r ( n ) ) and an+2(r(n)). 118 Chapter 8. Implementation & Error Discussion n r[n) M-(n),n, 15 an+i(r(n)) a n + 2 ( r (n)) 0 0.319264 5.5 x 10~3 +5.4 x 10~3 -7.7 x 10-3 1 1.489194 1.0 x 10"4 -1.0 x 10~4 +1.1 x 10"4 2 2.603209 6.3 x 10~7 +5.3 x 10~7 -3.4 x 10"7 3 3.655180 8.5 x 10~8 +8.4 x 10-8 -9.3 x lO" 8 4 4.340882 4.3 x 10-9 +4.2 x 10-9 -4.6 x 10"9 5 5.581000 1.2 x 10-1U -1.2 x 10"10 +1.2 x 10"1U 6 6.779506 2.7 x 10-1'2 +2.7 x 10"12 -2.5 x 10"12 7 7.879012 3.9 x 10-14 +3.6 x 10"14 -4.7 x 10"14 8 8.406094 6.9 x 10"15 +6.9 x 10~15 -7.1 x 10~15 9 9.656578 2.1 x 10"16 -2.0 x 10~16 +2.0 x 10"16 10 10.900511 6.i x n r 1 8 +6.1 x 10"18 -5.9 x 10"18 11 12.066012 1.1 x 10"19 - i . i x n r 1 9 +9.1 x 10~20 12 13.144565 5.2 x 10"21 -5.1 x 10-21 +5.6 x lO" 2 1 Table 8.6: Comparison of M r ( n ) n and an+\(r{n)) The data in Table 8.6 would lead one to conjecture that Mr(n),n ~ Wn+l(r(n))\ , which, if true, provides a much simpler method for bounding le , .^ , i n the right-half plane. 119 Chapter 9 Comparison of the Methods In this chapter we revisit the computation of T ( l + z) from Chapter 2 for the purpose of comparing Lanczos' method against those of Stir l ing and Spouge. For each of these methods T(20 + V7%) is computed wi th a relative error |e p | < 10~ 3 2 , and then formulas for T(l + z) are given for each wi th a uniform error bound of 10~ 3 2 in the entire right-half plane ReO) > 0. 9.1 Stirling's Series For this section recall the notation of Section 2.5: let r B2n(x) EN,n{z) = / ——dx , Jo 2n(z + N + x)2n denote the integral error term of equation (2.11), and denote by U^tn(z) the upper bound on Ej^^n(z) given by Theorem 2.3: ( 1 \ 2n+2 r> 1 \ £>2rc+2 cos (6/2) J (2n + 2) (2n + 1) (z + T V ) 2 ^ 1 ' 9 . 1 . 1 r ( 2 0 + 17i) with \ep\ < 1 0 - 3 2 Since z = 19 + 17i is fixed, for convenience write EN,U = EN^(\%-\-17i) and UNtn = UNj7l(19 + 17i). From equation (8.2), E ^ > N is related to the relative error e p v ia eEN,n _ 1 120 Chapter 9. Comparison of the Methods from which e p = 1 - e~EN'n We want |e p | < 10 3 2 which means |1 — e E N < n \ < 10 3 2 , which is guaranteed if \EN,N\ < \UN,n\ < 10" 3 2 / e « 3.6 x 1 0 " 3 3 . The pairs (JV, n) which achieve this bound are listed in Table 9.1. N n • \UN,TI\ < 18 10 3.5 x 10 3 3 10 11 2.9 x 10 3 3 4 12 3.3 x 10 3 3 .0 13 2.6 x 10 3 3 Table 9.1: Upper bound on |.Ejv,n(19 + lTz) j in Stirl ing Series Selecting (N, n) = (4,12), equation (2.11) yields log T(20 + I7i)l[(z + k=l (19 + 17^ + 4 + 1/2) log (19 + 17i + 4) - (19 + 17i + 4) + - log2vr 12 •B-2j ^ 2 j ( 2 j - l ) ( 1 9 + 17i + 4 ) 2 J - 1 ' from which T(20 + 17z) = -6.6530978807100357093202320786706 x 10 1 3 +i 1.3813486137818296429873066956513 x 10 1 4 . Comparing this value against a high precision value produced by Maple, we find a relative error of approximately 8.5 x 10~ 3 4 . The Bernoulli numbers in the sum range in size from B2 = 1.7 x 1 0 _ 1 to - B 2 4 = 8.7 x 10 4 . The corresponding terms of the sum range in size from 2.9 x 10~ 3 to 5.0 x 10~ 3 2 . 121 Chapter 9. Comparison of the Methods 9.1.2 T(l + z) with \ep\ < lO" 3 2 Uniformly The goal this time is a formula for T ( l + z) wi th relative error |e p | bounded uniformly by 1 0 - 3 2 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 \U^^n(z)\ becomes unbounded as z —> 0+ along the real axis. B y the same argument as in the previous section, again in this case the bound |C/jv,n(0)| < 3.6 x 10~ 3 3 is sufficient to ensure \ep\ < 10~ 3 2 . Computing |C/jv,n(0)I f ° r pairs (N,n) which achieve this bound results in Table 9.2. Selecting (N, n) = (17,17), for example, we N n | tW0)| < 19 15 3.5 x l f r 3 3 18 16 1.4 x 1 0 " 3 3 17 17 9.4 x 1 ( T 3 4 16 18 9.7 x l f r 3 4 15 19 1.7 x i r r 3 3 14 21 l . i x l f r 3 3 13 23 2.4 x l f r 3 3 12 28 2.2 x l f r 3 3 Table 9.2: Upper bound on \ENtn(0)\ in Stirl ing Series thus obtain the following formula log 17 r(i+ *)!!(* + *) /c=i (z + 17 + 1/2) log (z + 17) - {z + 17) + - log 2TT 17 +13 B 2j j^2j{2j-\){z + n ) ^ In this case, the Bernoulli numbers in the sum range in size from B% = 1.7 x 1 0 - 1 to = 4.3 x 10 1 1 , while the corresponding terms of the sum wi th z = 0 range in size from 4.9 x 10~ 3 to 9.5 x 10~ 3 3 . 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 y/2ir exp(—a) term which occurs in the leading fac-tor, and (27r) - 1 / 2 exp(a) which occurs in each of the coefficients (2.18). Spouge's formula may then be stated T(z + 1) = (z + a)z+1/2e-z where N = \a\ — 1, and N k=l ^z + k v ' (9.1) dk(a) = y^(-k + a)k-^e-k . (9.2) Under this simplification, the relative error in Spouge's formula remains the same: l e s ( a ' * ) l < (27r)^/ 'Re(* + a) " 9.2.1 T(20 + 17i) with \ep\ < 10 -32 For |e p | < 10 3 2 we require 'a ( 2 7 r ) a + V 2 R e ( l 9 + 17? + a ) < 10 -32 The left hand side of this expression is eventually decreasing as a func-t ion of a, and drops below 1 0 - 3 2 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 T(20 + I7i) = -6.6530978807100357093202320786706 x 10 1 3 - M 1.3813486137818296429873066956513 x 10 1 4 . In this instance, comparing this value against a high precision value produced by Maple, the relative error is approximately 2.6 x 10~ 4 4 , 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 <4(39) k dk(39) 0 +2.8947105233858572256681259383826xl0~ 1 7 20 - 1 . 4 6 1 2 0 4 3 5 9 6 4 4 5 0 4 4 3 8 9 3 4 0 5 4 0 6 6 3 5 9 2 X 1 0 - 1 1 +2 .2677611785616408448051623649413x10° 21 + 1.6856896402284750130501640213058xl0- 2 2 - 3 . 0 4 5 8 8 5 8 4 2 6 2 6 1 6 8 4 7 4 0 7 5 2 1 8 2 0 6 6 8 5 7 X 1 0 1 22 -1.5553541637142526870805344933539xl0- 3 3 +1.9357212181425501030368331204744xl0 2 23 +1.1302120506210549847428847358464 x 1 0 ~ 4 4 -7.7429949716371535150533416191734xl0 2 24 -6.3472060001627477262571023780708xl0- 6 5 +2.1876179361012802304828398825556xl0 3 25 +2.6919634338361067149994936318051 x l O " 7 6 - 4 . 6 4 3 8 5 3 7 2 7 7 5 4 8 5 6 0 4 1 7 6 2 3 3 6 0 4 0 5 7 7 5 X 1 0 3 26 -8.3801786194339100406682820681103xl0- 9 7 +7.6927388480783573540015702664876xl0 3 27 + 1.8481319798192137441154507903452 x l O - 1 0 8 -1.0195917484892786070501293142239xl0 4 28 -2.7610263895691596705265848094462 x l O ~ 1 2 9 +1.0999167217153749462567889811423xl0 4 29 +2.6382697777056113806964677703858 x 1 0 " 1 4 10 - 9 . 7 7 4 0 1 7 8 5 6 7 2 2 7 0 1 9 4 4 6 5 6 8 4 2 1 9 0 3 0 3 5 X 1 0 3 30 -1.4954804884394724082929843529081 x l O - 1 6 11 +7.2136820304845623606047714747868xl0 3 31 +4.5441806633655985113201441255734xl0- 1 9 12 -4.4462490659018128337839056551451x 10 3 32 ' - 6 . 4 2 8 9 9 4 7 4 5 1 7 2 2 9 5 2 2 6 3 5 3 4 4 2 4 1 5 9 0 5 3 X 1 0 - 2 2 13 +2.2961572174783000519663098692718xl0 3 33 +3.4516439287785093504202451110646xl0- 2 5 14 -9.9491746363493229520431441706026xl0 2 34 -5.1380351805226988204329822384216 x 1 0 " 2 9 15 +3.6160749560761185688455782354331xl0 2 35 + 1 . 2 6 0 6 6 0 7 4 6 1 7 8 7 9 7 6 9 4 3 2 3 4 8 1 1 1 8 5 4 4 6 X 1 0 " 3 3 16 -1.1004495488313215123003718194264xl0 2 36 -1.9452281496812994762288072286241 x l O - 3 9 17 + 2 . 7 9 4 7 8 5 2 8 1 6 1 9 5 5 7 8 9 9 4 2 0 4 7 3 1 1 9 5 1 5 4 X 1 0 1 37 +2.2292761113182594463394149297371 x l O - 4 7 18 -5 .8947917688963354031891858182660x10° 38 -2.2807244297698558308437038776471x 1 0 ~ 6 ° • 19 +1 .0259323949497584239156793246602x10° Table 9.3: Coefficients of Spouge's series, a = 39 9.2.2 T( l + z) with \ep\ < 10~32 Uniformly In the case of Spouge's method wi th a uniform bound the procedure is essentially the same. The bound on the relative error is greatest when Re(^) = 0, so for \ep\ < 10~ 3 2 we require 1 . l n - 3 2 - < 10" ( 2 7 r ) a + l / 2 a Aga in 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 e^n, and determine the least n which yields the prescribed error bound. For |e p | < K T 3 2 we find using Table C . l that n = 21, r = 22.618910 gives the empirical bound \ertTl(it/(l — t))\ < 2 x 1 0 - 3 4 for t G [0,1). Thus, by equation (8.6), W < Q " 2 x l O -< 2.2 x 1 0 " 3 4 . Using the formulation (6.14), we find r(, + i ) - 2 , / « ' * + r<21> + W * H / ' 7T v e 2 1 J t 1 z + k where the coefficients dk are given in Table 9.4. Evaluating this expression at z = 19 + 17?; we find to 32 digits T(20 + 17i) = -6.6530978807100357093202320786706 x 10 1 3 + i 1.3813486137818296429873066956513 x 10 1 4 . The relative error in this calculation is less than 1.5 x 10~ 4 2 . 9.4 Discussion A l l 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 4 0 +2.0240434640140357514731512432760 x 1 0 " 1 0 1 +1.5333183020199267370932516012553 x 10° 2 -1.1640274608858812982567477805332 x 10 1 3 +4.0053698000222503376927701573076 x 10 1 4 -8.2667863469173479039227422723581 x 10 1 5 +1.1414465885256804336106748692495 x 10 2 6 -1.1135645608449754488425056563075 x 10 2 7 +7.9037451549298877731413453151252 x 10 1 8 -4.1415428804507353801947558814560 x 10 1 9 +1.6094742170165161102085734210327 x 10 1 10 -4.6223809979028638614212851576524 x 10° 11 +9.7030884294357827423006360746167 x 10" 1 12 -1.4607332380456449418243363858893 x 10" 1 13 +1.5330325530769204955496334450658 x 10" 2 14 -1.0773862404547660506042948153734 x 10~ 3 15 +4.7911128916072940196391032755132 x 10" 5 16 -1.2437781042887028450811158692678 x 10" 6 17 +1.6751019107496606112103160490729 x 10~ 8 18 -9.7674656970897286097939311684868 x H T 1 1 19 +1.8326577220560509759575892664132 x 10~ 1 3 20 -6.4508377189118502115673823719605 x 1 0 " 1 7 21 +1.3382662604773700632782310392171 x 1 0 " 2 1 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. A l l 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 T ( l + z) wi th a uniformly bounded error, the Lanczos method appears to be the clear winner. In the case \ep\ < 10~ 3 2 considered here, 21 divisions and 21 additions were required in the 126 Chapter 9. Comparison of the Methods sum. B y 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+V wi th increasing \z\. Spouge's error decreases wi th increas-ing \z\, but only like |Re(z + The Lanczos error also decrease to zero as \z\ —> oo, by virtue of the choice of r as the largest zero of e^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 l o g r ( l + z) as in Stirling's series, and then take the exponential of the result. A s 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. A s was noted in Section 2.5.2, it is possible to achieve a relative error of order 10~ 3 4 in the calculation of T(7 + 13z) using (N,n) =-(0,38), but the largest Bernoulli number required wi l l be of order 10 5 0 . Roundoff can also be a problem in the summation of the series themselves, even wi th accurately computed coefficients. For example, the coefficients in Table 9.3 span 63 orders of magnitude, and many significant digits of some individual terms wi l l 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 [—TT/2, TT/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 Limi t Formula. Finally, an improvement to Stirling's formula is noted based on the existence of zeros of the error function e£°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, in (3.6), r (^ + l / 2 ) ( ^ + r + l /2 ) - ( z + 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 128 Chapter 10. Consequences of Lanczos' Paper where g(6) € L2 [-TT/2, TT/2] is even. Note that g is very general, the only requirement being that it be square summable wi th respect to Lebesgue measure. Then the Fourier coefficients 2 r / 2 ak = - g{9) cos (2fc0) dO . are defined, and we can associate with g its Fourier series oo 9(0) ~ 5^ ' °fe c o s ( 2 k e ) • The ~ symbol is used since the series converges to g in L2[—ir/2,TT/2], though not necessarily pointwise. Nonetheless, since cos 2 2 6 £ L2[—n/2, n/2], it follows from Parseval's theorem and identity (3.9) of Lemma 3.2 that F{z) = / cos 2 2 6g{&) d9 J-ir/2 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 i n the proof of Theorem 5.1 shows that Hk(z) = O ( A T 1 / 2 - " 5 ) on K for some S > 0. B y the Schwarz inequality, (±\at\mz)X<c(±\aTA (E^U \fe=AT / \k=N / \k=N J 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 AT 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^(9) is continuous wi th g(-k\—7r/2) = g^k\ix/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 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) wi l l converge absolutely and uniformly on K, and thus defines an analytic function in the region {z e C | Re(z) > - m / 2 and z ^ - 1 , - 2 , - 3 , . . . } . Since the r(2+l/2)/T(z+l) 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-t ion 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 the series in (10.1) after the nth term, the resulting error er<n(z) = 2^feLn+i akHk(z) has the property that in the half plane Q, = Re(z) > 0, a constant. Invoking Lemma 8.2, the maximum of | e r , n ( ,2 ) | in Q occurs on the line Re(z) = 0, and the maxima'can be located and estimated by examining the first few terms of series \erjTl(it/(l — t))\, 0 < t < 1. A s 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-t ion 6.1 carry over word for word, wi th the Horner method probably the more practical. If g{9) 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-oo n - Ylak 130 Chapter 10. Consequences of Lanczos' Paper als, Tk(cos0) = cos (kO). Thus 2 fn/2 ak = - g{6) cos {2k9)d0 71 J-n/2 2 pn/2 k = Z 9(9) J2C^ COS2J DDD = lT,C^ I ' COS 2'00(0) d0 = lJlCmkF{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 2 k a>k = - Y C2j,2kF(j) j=0 and form the formal sum oo g(9) = ^ ak COS (2kd) . k=0 If this series is in L2[—n/2,n/2] , that is, if J2T=oal < 0 0 > then it defines g{9) as the inverse transform of F(z). An Example The following example is found in [13, p.45]. The Bessel function of the first k ind of order z not necessarily an integer can be expressed i (x/2)z r/2 J ^ x ) = -7= TV i i ^ / c o s 2 ^ c o s ( x s i n 0 ) d 0 . V7rl (z + 1/2) J_n/2 131 Chapter 10. Consequences of Lanczos' Paper Lett ing F(z,x) = V^ Jz(x)T(z + l/2){x/2)~z , we see that F(z, x) is the Lanczos transform of cos (xsrnO). Hence -T(z + 1/2) 00 F(z,x) = ^ { * z ' Y - a k { x ) H k [ z ) k=0 so that = T^TTTE' ' ^ (x )H k (z ) 1 ^ + l) k = 0 The coefficients ak(x) are given by 2 k ak(x) = - V] C2j,2kF(j,; 2 fc -7= X > ; , 2 f c J , ( z T C + l/2)(x/2)~i j=0 O n the other hand, wi th the transformation 0 = 9 —TT/2, the Fourier coefficients become 2 / ' 7 r / 2 «fe (^) = — / cos (x sin 9) cos (2A;#) d# T J-K/2 2 r = - / cos (z sin(</> + TT/2)) cos (2fc(^.+ TT/2)) d<£ Ti" Jo 2 / > 7 r = (—l)k— / cos (x cos (p) cos (2kcf>) dcf> Tr Jo = 2J2fc(^ 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: 1 k J2k(x) V5F Y2CwJj{x)rti + l/2)(x/2)-i . More importantly, 2{x/2)z T(z + 1) oo YJ2k(x)Hk(z) . fc=0 It is interesting in this case to examine what the properties of g(9) = cos (x sin 9) say about the rate of convergence of the series. In this case, ^r ( f c ) (—7r /2 ) = g( f e)(7r/2), for al l k > 0, so that for fixed x, ak(x) = 2J2k{x) = O (k~n) for any n > 0 as k —> oo. Furthermore, again for fixed x, Jz(x) is entire as a function of z. The following non-trivial identities are consequences of the Lanczos Limi t Formula 7.1: L e m m a 10.1. Suppose n > 1 is an integer. Then 10.2 A Combinatorial Identity (ix) P r o o f of L e m m a 10.1: From 7.1 1 oo T(z + 1) = 2 l im rz Setting z = n terminates the series giving n 133 Chapter 10. Consequences of Lanczos' Paper Since rn —>• oo while nl remains constant, it follows that the term in square brackets must tend to zero as r —> oo. This gives (i). For (ii), again setting z = n terminates the series which may be writ ten 1 I N p n / - i \ k r - k 2 / r n\n\ . 0 2 _ r Z ^ f c = l l J-J e (n-k)\{n+k)\ nl = 2 l im = — — . i—>oo r n Apply ing L 'Hopi ta l ' s rule n times then gives En i -t\kr-k2/r win! L2 fc=iv ^ e fn-fc)!(n+/c)!ft' n! = 2 l im _ n r - n + l En f 1 \k„-k2/r n\n\ / 4 fc=ll e (n-fc)!(n+fc)!ft' £ 2 l im (n-fc)!(n+fc) r-»oo n(n — 1 ) ? — n + 2 E n { l\kr-k2/r n\n\ ]2n = 2 l im {n-k)\{n+k)\K r—>oo n ( - l ) n n! ^ l j ( n - f c ) ! ( n + A:)! 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 T(z + 1) = e - ^ 2 + 1 / 2 ( 2 7 r ) 1 / 2 e ^ m ' (10.2) for some 0 < 6 < 1, wi th resulting relative error exp (#/(12z)) — 1 « #/(12z) which decreases wi th increasing z. However, without some manipulation, the formula is not very effective for small z. Spouge in [27] derives the more visually pleasing approximation T(z + 1) « (z + 1 /2 ) 2 + 1 /2 e - (^+i /2)^ ( 1 0 3 ) with bounded relative error ^ J v^Fi TT FF + l/2 ' 134 Chapter 10. Consequences of Lanczos' Paper which is at most 0.107 at o = Re(z) = 0. A s 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 + l / 2 i s introduced to simplify the resulting formulas. The work in Section 8.5.4 furnishes an improvement to (10.2), albeit wi th the introduction of a slightly complicated constant. Let t ing r(0) denote the largest zero of e^Q, and setting a = r(0) + 1/2, we find T(z + l) = (z + a ) ( * + 1 / 2 ) e - ( * + a ) 7> H er(0),0(2) since 0 = 0 implies ao(r(0))/2 = 1. This approximation is valid in the entire right half plane Re(z) > 0, is asymptotic to T(z + 1) as \z\ —>• oo in this region, and has an empirically determined relative error |er(0),o(*)| < 0.006. The constant a = 0.819264 appears complicated, but it turns out that the solutions of e£°n = 0 can be found explicitly in terms of Lambert W functions in the n = 0 wi l l now be shown. L e m m a 10.2. Let r(0) be the largest zero o /e^° 0 ; and a = r(0) + 1/2. Thena = -W(-l/n)/2. P r o o f o f L e m m a 10.2: Recall that 2 ~ ^ ( r + 1/2) y/2na 135 Chapter 10. Consequences of Lanczos' Paper The real zeros of e£°0 are the solutions of 1 — ao/2 = 0, so V27ra 27ra = e 2 a - 2 a e " 2 a = - - . 7T That is, —2a = W(—l/n) whose only real roots are a=-W0(-l/7r)/2 = 0.276914 and a=-W-1(-l/ir)/2 = 0.819264 . • The value of a = — W^\(—l/n)/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 T{z + 1) = (z + a{z)Yz+1^e-^z+a^V2^ ? Recalling the definition of Fr(z) from Section 8.2 Fr(z) - r ( z + l)(z + r + l / 2 ) - ( 2 + 1 ^ e 2 + r + 1 / 2 ( 2 7 r ) - 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 oi.Fr(z) w i l l be needed (see [27]): 136 Chapter 10. Consequences of Lanczos' Paper Lemma 10.3. F0(z) strictly increases to 1 as z —> oo, and Fo(0) = 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) E (—1/2, oo) x [0, oo). Secondly, observe that l im Fr(z) — 1 = l im Fr(z) — 1 = oo . r—1—1/2 r^oo Thirdly, by Lemma 10.3, FQ(Z) — 1 ] 0 as z —» oo, so that,Fo(z) — 1 < 0. Thus for z > 0 fixed, /^.(z) has at least two zeros as a function of r . A n explicit formula for a(z) can be found wi th the help of the Lam-bert W function. Suppose that z is real and positive. Then so that T(z + l) = (z + a ( z ) ) ( z + 1 / 2 ) e - ( 2 + a ( z ) ) V ^ T T , T ( z + 1 ) = (z + a(z))l'+1'%-l'+aW . 2TT Raising both side to the power 1/(3+1/2) and then multiplying through by - 1 / ( 3 + 1/2) gives T(z + 1) z + l/2[ ^/2rr whence z + a(z) ' * + l/2 1 / ( 2 + 1 / 2 ) 3 + 0 ( 3 ) z + 1/2 r(* + i) and a(z) = - z - ( z + 1/2) VFfe + 1/2 - 1 3 + 1/2 2rr 1 / ( 2 + 1 / 2 ) ^ r(3 + i) 2vr l / ( Z + l / 2 ) ^ (10.5) 137 Chapter 10. Consequences of Lanczos' Paper For real z the branches W _ i and WQ 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 V itself, and the question of whether these equations gives rise to more efficient T approximation schemes was not studied further. Equation (10.5) does, however, give a picture of what a(—l,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/(l — z), 0 < z < 1. It is interesting to observe how little these function appear to vary over the entire real line, in particular a(—l,z) which is approximately 0.82 at z = 0 and decreases to about 0.79 as z oo. ' •^- i (^ ) j • >a(— 1,2 /(I-*) ) j 1 j a(a.z/( | Figure 10.1: Pa th of o(0, jzrz), a(-l, 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: W i t h respect to Lanczos' paper itself: 1. Wha t is the precise asymptotic bound on the coefficients afe(r) of the main series Sr(z) as k —> oo? 2. C a n 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 L imi t Formula (7.1) be extended to Re(z) < 0? W i t h respect to the error bounding methods of Chapter 8: 1. Is there an analytical method to bound the relative error function k r , n ( ^ / ( l — t))\ for 0 < t < 1, and hence | e r > n ( £ ) | uniformly in the right-half plane? More simply, in the n = 0 case, is there an analytic bound for |e r ) 0 ( i t ) | = \Fr(it) — 1|? 2. Wha t is the precise relationship between the largest zero of the relative error at infinity, e^n, and the maximum of \erjTl(it/(1 — i ) ) | on 0 < t < 1? 139 Chapter 11. Conclusion and Future Work 3. How many zeros does e£°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, wi th respect to a deterministic algorithm: 1. Is there a deterministic algorithm which, given z and e > 0, pre-scribes n and r so that T(z + 1) is computed wi th relative error at most e? 2. Wha t 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 a/c(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 —> oo. Lanczos' state-ments suggest, and numerical evidence supports, the slightly faster rate ak(r) = 0(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~l) as k —• oo, 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) = T(z + l)(z + r + l / 2 ) - z - 1 / 2 e 2 + r + 1 / 2 ( 2 7 r ) - 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 unti l 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 al l z € C away from the negative integers. His claim appears to be based on the asymptotic behaviour of the coefficients ak(r) as r —> oo. Indeed, convergence of the series Sr(z) extends further and further left in the complex plane as r increases, but just how oo l im (z + r + l / 2 ) z + 1 / 2 e ^ z + r + 1 ^ V' ak(r)Hk{z) r—>oo ^—' fc=0 CO = 2 l im rz Y\-l)ke-k2/rHk(z) r—>co ^—' fc=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 er,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 er<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 al l cases examined, setting r equal the largest zero of e^°n results in uniform error bounds much improved on Lanczos' own estimates. W h y 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 er,n(z) as a function of r wi th 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 er,n(z) = 0, and consequently the approximation to F(z + 1) becomes exact. Empir ica l 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^(z), —0.5 < r < 3, —0.5 < z < 3 which shows several zeros of er^(z) for each value of z. For Re(z) > —1/2 and n > 0 a fixed integer, if r(z) 3 3 Figure 11.1: e r, n(^), - 0 . 5 < r < 3, - 0 . 5 < z < 3 is defined to be a zero of ertn(z), equation (10.5) and Figure 11.1 lead one to conjecture that r(z) is a multivalued complex analytic function wi th number of branches equal to the number of zeros of e^n. The data from Tables C . l and C.2 suggests that neither the number nor the distribution of e£°n zeros follows a simple law. Even more per-plexing is the striking similarity between the functions and an+i(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 al l 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 wi th Re (2) > 0 and e > 0 as input, prescribes optimal n and r and computes T(l-\-z) wi th relative error at most e. A s noted, the Lanczos method lacks a simple functional relationship between the relative error er<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 wi th a uniform error bound, the problem is st i l l difficult. In the uniform case, however, the empirical evidence strongly sug-gests, a simple deterministic algorithm. This evidence is the data of Tables C . l 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 T(l + z) wi th uniform relative error e > 0 requires O (— loge) terms of the sum, as given by the bound (2.19). In the Lanczos plot of pairs (n, — l o g M r ( 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 1 o g M r ( n ) , n . Thus given e > 0, one should choose n = \—2 — 0.3 log e]. The corre-sponding r should then be chosen to be the largest zero of e^ , which is about size n. This choice of r « n is consistent wi th 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 Bibliography T . Apostol . Introduction to Analytic Number Theory. Springer, New York, 1998. E . A r t i n . The Gamma Function. Holt , Rinehard and Winston, New York, 1964. J .P. Boyd. Chebyshev and Fourier Spectral Methods. Dover Pub-lications, Inc., New York, 2000. F . Cajori . A History of Mathematical Notations. The Open Court Publishing Company, Chicago, 1929. R . M . Corless, G . H . Gonnet, D . E . G . Hare, D . J . Jeffrey, and D . E . Knu th . O n the Lambert W function. Adv. Comput. Math., 5(4):329-359, 1996. P. J . Davis. Leonhard Euler's integral: A historical profile of the gamma function. The American Mathematical Monthly, 66:849-869, 1959. Serge Dubuc. A n approximation of the gamma function. J. Math. Anal. Appi, 146(2):461-468, 1990. H . M . Edwards. Riemann's Zeta Function. Academic Press, 1974. Pau l Godfrey. A note on the computation of the conver-gent lanczos complex gamma approximation. Web published at h t t p : / / m y . f i t . e d u / ~ g a b d o / p a u l b i o . h t m l , 2001. D . Gronau. W h y is the gamma function so as it is? TEACHING MATHEMATICS AND COMPUTER SCIENCE, 1:43-53, 2003. C . Lanczos. Trigonometric interpolation of empirical and analyt-ical functions. Journal of Mathematics and Physics, 17:123-199, 1938. C . Lanczos. Applied Analysis. Prentice Hal l , New Jersey, 1956. C . Lanczos. Linear Differential Operators. D . V a n Nostrand Com-pany, L t d . , 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. Nor th Carolina State University, College of Physical and Mathematical Sciences, Raleigh, N C , 1998. W i t h a preface by W i l l i a m R. Davis, a foreword by George M a r x and a biograph-ical essay by Barbara Gellai , Translated by Jozsef Illy, Laurent Choquet, Don Ridgeway and Judi th Kontsag Mesko. [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 & H a l l / C R C , Boca Raton, U S A , 2003. [19] P. Mikusinski and D . Taylor. An Introduction to Multivariable Analysis from Vector to Manifold. Birkhauser, Boston, U S A , 2002. [20]' Edward W . Ng . A comparison of computational methods and algorithms for the complex gamma function. ACM Trans. Math. Software, l ( l ) :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, U K , 1992. [22] T . J . R i v l i n . Chebyshev Polynomials from Approximation Theory to Algebra and Number Theory. John Wi ley & Sons, Inc., New York, 1990. [23] W . Rudin . Principles of Mathematical Analysis. M c G r a w - H i l l , Inc., New York, 1976. [24] W . Rudin . Real and Complex Analysis. M c G r a w - H i l l , 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. Computat ion of the gamma, digamma, and tr igamma functions. SI AM 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 Wi ley & 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 . Whit taker and G . N . Watson. A Course of Modern Analy-sis, Fourth Edition. Cambridge University Press, Cambridge, U K , 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 wi th (3.1.1). From poo T(z + 1) = az+1 / tze~atdt , Re(a) > 0 , Jo make the substitution a = 1 + pz , where p > 0. This gives poo T{z + 1) = (1 + pz)z+1 / [te^Ye-'dt . Jo Next introduce the factor ( e p ) - 2 poo T{z + 1) = (1 + pz)z+1{ep)-z / (pte^yJdt , Jo and make the replacement v = el~pt to get T(z + 1) = (1/p + z ) * + 1 e - z - 1 / " / [v(l ~ log v) ] 2 u 1 / " - 1 ^ . Jo Finally, let r = 1/p — 1 to arrive at (z + r + i ) * + i e - ( * + H - i ) f _ logv)]* , Jo 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 wi l l be concerned only wi th this first kind, which is the only polynomial solution for integer n of the differential equation (i-^s-«i+»^=o. The defining relationship for the {Tn(x)} is Tn(x) = cos (nO) where x — cos# , and where x ranges from — 1 to 1 as 0 ranges from 7r to 0. The resulting graph of Tn(x) oscillates between —1 to 1 and looks like a distorted version of'cos (n6>) . The leading coefficient of Tn(x) is 2 n _ 1 and the function is even or odd according as n is even or odd. The first few Chebyshev polynomials of the first k ind are: T0(x) = 1 Ti(x) = x 149 Appendix B. A Primer on Chebyshev Polynomials T2(x) = 2x2 - 1 T3(x) = 4 z 3 -3x T 4 (x ) = 8 x 4 - 8 x 2 + l T5(x) = 1 6 x 5 - 2 0 x 3 + 5 x . The family {Tn(x)}^=0 is orthogonal with respect to the inner prod-uct { f ' 9 ) = J j { x ) 9 { x ) 7 i = h ( R 1 ) and forms a basis for C2[— 1,1] wi th 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 n=l where convergence of the series is uniform. The coefficients Cn(r) are given by cn(r) = - f f(x)Tn(x)^L= . (B.2) TT J_i V I - X2 A s an operational detail, take note that the family {Tn(x)}'^>=0 is not an orthonormal set, but only orthogonal. The proper normalization is { v^r°' ^ T n { 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 0 which can then be written as a Fourier cosine series. A s 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(VT^) = (-l)nT2n(x) (B.3) 150 Appendix B. A Primer on Chebyshev Polynomials and the trigonometric form T 2 n (s in#) = ( - l ) n c o s ( 2 n 0 ) . (B.4) Finally, the following important property of Tn(x) should be noted: Of al l monic polynomials of degree n o n [-1,1], the one wi th 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 . l and C.2 summarize the results of numerical investigations carried out on the relative error functions e r > n (z) and e ^ . For each n = 0 , . . . , 60, the number of zeros of e^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 er>n(it), 0 < t < oo based on the theory of Sec-t ion 8.5. For r equal the largest zero of e£°n, the maximum of was estimated using the first 5 and then 15 terms of the sum. These columns are labeled Mr(n^n<5 and M r ( n ) j n ) i 5 , respectively. In addition, the location of the maximum of ak{r)Hk{tt)\ was determined; this data is in the last column labeled t m a x . 152 Appendix C. Empirical Error Data n No. Zeros Smallest Zero Largest Zero r(n) •Wr(n),n,5 -Mr(n),n,15 ^max 0 2 -.223086 .319264 5.9 x 10~ 3 5.5 x 10" 3 .88 1 4 -.173495 1.489194 1.0 x 10~ 4 1.0 x 10" 4 2.13 2 6 - .151082 2.603209 6.4 x 10~ 7 6.3 x 10" 7 4.27 3 8 - .137917 3.655180 8.5 x 10~ 8 8.5 x 10" 8 4.05 4 8 - .129067 4.340882 4.3 x 10" 9 4.3 x 1 0 - 9 5.05 ' 5 10 -.122605 5.581000 1.2 x 1 0 - 1 0 1.2 x 10~ 1 U 6.34 6 12 - .117620 6.779506 2.7 x 10~ 1 2 2.7 x 1 0 - 1 2 7.88 7 14 -.113619 7.879012 3.9 x 1 0 - 1 4 3.9 x 1 0 " 1 4 6.72 8 14 -.110313 8.406094 6.9 x 1 0 - 1 5 6.9 x 1 0 " 1 5 9.23 9 16 -.107519 9.656578 2.1 x 10~ 1 6 2.1 x 1 0 " 1 6 10.51 10 18 -.105114 10.900511 6.1 x 10~ 1 8 6.1 x H T 1 8 11.83 11 20 -.103013 12.066012 1.1 x 10~ 1 9 1.1 x 1 0 " 1 9 14.30 12 22 - .101157 13.144565 5.2 x 10~ 2 1 5.2 x 1 0 - 2 1 12.38 13 22 - .099499 13.726821 4.0 x 10~ 2 2 4.0 x 10~ 2 2 14.69 14 24 -.098005 14.977863 1.2 x 1 0 - 2 3 1.2 x 1 0 " 2 3 15.98 15 26 - .096650 16.209805 3.6 x 1 0 - 2 5 3.6 x 10~ 2 5 17.45 16 28 -.095412 17.345444 3.1 x r o - 2 7 3.1 x 10~ 2 7 23.47 17 30 -.094275 18.399283 5.0 x 1 0 - 2 8 5.0 x 1 0 " 2 8 18.11 18 30 -.093226 19.048512 2.5 x 1 0 " 2 9 2.5 x 10~ 2 9 20.16 19 32 -.092252 20.298892 7.8 x 10~ 3 1 7.8 x 1 0 - 3 1 21.45 20 34 -.091346 21.508926 2.1 x 10~ 3 2 2.1 x 10~ 3 2 23.40 21 36 -.090499 22.618910 1.8 x 1 0 - 3 4 1.8 x 1 0 - 3 4 17.14 22 36 -.089704 23.118012 5.2 x 10~ 3 5 5.2 x 1 0 " 3 5 24.34 23 38 -.088958 24.370498 1.7 x 10~ 3 6 1.7 x 1 0 - 3 6 25.62 24 40 -.088253 25.617904 5.2 x 1 0 - 3 8 5.2 x 1 0 - 3 8 26.97 25 42 -.087588 26.798597 1.1 x 1 0 - 3 9 1.1 x 1 0 " 3 9 30.13 26 44 - .086957 27.886311 3.6 x 1 0 - 4 1 3.6 x 1 0 " 4 1 24.97 27 44 -.086358 28.440357 3.5 x 1 0 - 4 2 3.5 x 1 0 " 4 2 29.80 28 46 -.085789 29.692534 1.1 x 1 0 - 4 3 1.1 x 1 0 " 4 3 31.09 29 48 -.085246 30.931341 3.4 x 1 0 - 4 5 3.4 x 1 0 " 4 5 32.63 30 50 - .084727 32.080670 4.4 x 1 0 " 4 7 4.4 x I Q " 4 7 39.41 Table C . l : Error Data, n = 0 , . . . . 30 153 Appendix C. Empirical Error Data n No. Zeros Smallest Zero Largest Zero r(n) M-(n),n,5 •Wr(n),n,15 m^ax 31 52 - .084232 33.145772 4.2 x 1 0 " 4 8 4.2 x 1 0 " 4 8 31.63 32 52 -.083757 33.762726 2.4 x l O " 4 9 2.4 x 10~ 4 9 35.27 33 54 - .083302 35.014250 7.7 x 10~ 5 1 7.7 x 1 0 " 5 1 36.57 34 56 - .082865 36.235367 2.2 x 10~ 5 2 •2.2 x 10~ 5 2 38.65 35 58 - .082445 37.356480 7.5 x 10~ 5 5 7.5 x 10~ 5 5 34.71 36 60 -.082041 38.385241 3.8 x 1 0 " 5 5 3.8 x 1 0 - 5 5 38.30 37 60 -.081651 39.085095 1.7 x l O " 5 6 1.7 x 1 0 " 5 6 40.74 38 62 - .081275 40.334630 5.3 x 10~ 5 8 5.3 x 1 0 " 5 8 42.09 39 64 - .080912 41.529155 1.3 x l O " 5 9 1.3 x 10~ 5 9 45.44 40 66 -.080562 42.626437 2.7 x 10~ 6 1 2.7 x 10~ 6 1 36.11 41 66 -.080223 43.154830 3.6 x 1 0 " e 2 3.6 x l O " 6 2 44.92 42 68 -.079894 44.407411 1.2 x 1 0 " 6 3 1.2 x 10~ 6 3 46.21, 43 70 -.079576 45.651117 3.7 x 10~ 6 5 3.7 x l O " 6 5 47.74 44 72 -.079268 46.814382 6.2 x 1 0 " 6 7 6.2 x 10~ 6 7 54.32 45 74 -.078969 47.889652 3.8 x l O " 6 8 3.8 x 1 0 " 6 8 44.46 46 74 -.078679 48.477371 2.6 x 1 0 " 6 9 2.6 x l O " 6 9 50.39 47 76 -.078396 49.729491 8.2 x 1 0 " 7 1 8.2 x 1 0 " 7 1 51.68 48 78 -.078122 50.959691 2.5 x 10~ 7 2 2.5 x l O " 7 2 53.73 49 80 -.077855 52.092791 1.7 x 1 0 " 7 4 1.7 x l O " 7 4 71.81 50 82 -.077596 53.141340 3.8 x 1 0 " 7 5 3.8 x lO"" 7 5 51.74 51 82 -.077343 53.799879 1.8 x l O " 7 6 1.8 x 1 0 " 7 6 55.85 52 84 -.077096 55.050733 5.8 x l O " 7 8 5.8 x l O " 7 8 57.19 53 86 -.076856 56.257932 1.5 x l O " 7 9 1.5 x l O " 7 9 60.47 54 88 - .076622 57.365268 1.8 x 1 0 " 8 1 1.8 x 1 0 " 8 1 44.65 55 88 -.076393 57.869538 4.0 x l O " 8 2 4.0 x 1 0 " 8 2 60.03 56 90 -.076170 59.122331 1.3 x l O " 8 3 1.3 x 10~ 8 3 61.32 57 92 -.075952 60.369399 4.1 x 10~ 8 5 4.1 x 1 0 " 8 5 62.81 58 94 - .075739 61.546699 8.5 x l O " 8 7 8.5 x l O " 8 7 68.92 59 96 -.075531 62.631604 3.3 x l O " 8 8 3.3 x 1 0 " 8 8 56.29 60 96 - .075327 63.192152 2.9 x 10~ 8 9 2.9 x l O " 8 9 65.50 Table C.2: Error Data, n = 3 1 , . . . , 60 154 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080001/manifest

Comment

Related Items