COMPUTING LAPLACE TRANSFORMS FOR NUMERICALINVERSION VIA CONTINUED FRACTIONSJoseph Abate and Ward Whitt900 Hammond Road, Ridgewood, NJ 07450-2908AT&T Labs, Shannon Laboratory, Florham Park, NJ 07932-0971; wowresearch.att.comINFORMS J. Computing 11 (1999) 394405AbstractIt is often possible to effectively calculate proba-bility density functions (pdfs) and cumulative distributionfunctions (cdfs ) by numerically inverting Laplace trans-forms. However, to do so it is necessary to compute theLaplace transform values. Unfortunately, convenient ex-plicit expressions for required transforms are often unavail-able for component pdfs in a probability model. In thatevent, we show that it issometimes possible to find continued-fractionrepresentationsfor required Laplacetransformsthatcan serve as a basis for computing the transform valuesneeded in the inversion algorithm. This property is verylikely to prevail for completely monotone pdfs, becausetheir Laplace transforms have special continued fractionscalled S fractions, which have desirable convergence proper-ties. We illustrate the approach by considering applicationsto compute first-passage-time cdfs in birth-and-death pro-cesses and various cdfs with non-exponential tails, whichcan be used to model service-time cdfs in queueing models.Included among these cdfs is the Pareto cdf.Keywordscomputationalprobability, numerical transforminversion, continued fractions, Laplace transforms, S frac-tions, complete monotonicity, Pade approximants, cumula-tive distribution function, birth-and-death process, Paretodistribution.Many probability density functions (pdfs) and cumu-lative distribution functions (cdfs) of interest in queue-ing models and other probability models arising in opera-tions research can be effectively computed by numericallyinverting Laplace transforms; see Abate, Choudhury andWhitt 1, Abate and Whitt 4, 5 and references therein.The biggest challenge in this approach, when there is achallenge, is usually computing the required Laplace trans-form values, because convenient closed-formexpressionsforLaplace transforms often are not available. In this paperwe point out that continued fractions can sometimes serveas a basis for effectively computing the required Laplacetransform values needed in the inversion algorithms.A simple motivating example is the steady-state waiting-time pdf in the M/G/1 queue. The classical Pollaczek-Khintchine (transform) formula gives the Laplace trans-form of the steady-state waiting-time pdf in terms of theLaplace transform of the service-time pdf. Thus we cancompute the waiting-time transformvalues in orderto com-pute the waiting-time pdf or cdf by numerical inversionwhenever we can compute the service-time transform val-ues. A possible difficulty, however, is that we might wantto consider service-time pdfs for which convenient explicitexpressions for the Laplace transform are unavailable. In-deed, this difficulty often arises when we consider distribu-Subject classifications: Probability distributions: calculation by trans-form inversion. Queues, algorithms: Laplace transform inversion. math-ematics, functions: Laplace transforms and continued fractionstions which have non-exponential tails, e.g., which cannotbe represented as phase-type distributions. The presentpaper providesa way to address this problem: Under favor-ablecircumstances, wemaybe abletoconstructacontinued-fraction representation of the service-time Laplace trans-form that enables us to compute the service-time Laplacetransform values, which in turn enables us to compute thewaiting-time Laplace transform values needed to performthe desired numerical inversion. A specific example coveredby this approach is the Pareto pdf.For background on continued fractions and their use fornumerical computation, see Baker and Graves-Morris 12,Bender and Orszag 13, Chapter 12 of Henrici 26, Jonesand Thron 28, Section 5.2 of Press, Flannery, TeukolskyandVetterling 32andWall 35. Applicationsof continuedfractions in statistics and applied probability are describedin Bowman and Shenton 15 and Bordes and Roehner 14.More recently, Guillemin and Pinchon 20, 21, 22, 23have used continued fractions to analytically derive impor-tant properties of queueing models. A summary of thatwork is contained in Dupuis and Guillemin 16. However,continued fractions evidently have not been suggested pre-viously as a way to numerically compute transform valuesin order to perform numerical transform inversion.The use of continued fractions is an alternative to com-putation of Laplace transforms via infinite-series represen-tations, which we recently discussed in Abate and Whitt9. We make an explicit numerical comparison to showthat continued fractions can be far superior in some cir-cumstances, even when the series converges geometrically.(See Section 6.)Here is how the rest of this paper is organized: In Sec-tion 1 we briefly define continued fractions and s