We present a simple, fast, and robust algorithm for numerically computing the first N moments (arbitrary N) of a nonnegative probability distribution from its Laplace–Stieltjes transform (continuous-mixed case) or z-transform (discrete case). The algorithm is based on numerically inverting an adaptively modified moment generating function. It only requires computation of the transform as several complex values of its argument. We also show that the high-order moments may be used in detecting the presence of exponential or geometric tails in distributions, and in case they are present the two parameters characterizing such a tail may be accurately computed. Several numerical examples of interest in the queueing literature illustrate the use of the algorithm. They include commonly used distributions as well as the waiting time, queue length, and busy period in queues with Poisson or more general Markovian arrival processes. Priority queues are also considered.
