  Stable distributions are an important class of infinitely-divisible
probability distributions, of which two special cases are the Cauchy
distribution and the normal distribution. Aside from a few special cases, the
density function for stable distributions has no known analytic form, and is
expressible only through the variate's characteristic function or other
integral forms. In this paper we present numerical schemes for evaluating the
density function for stable distributions, its gradient, and distribution
function in various parameter regimes of interest, some of which had no
pre-existing efficient method for their computation. The novel evaluation
schemes consist of optimized generalized Gaussian quadrature rules for integral
representations of the density function, complemented by various asymptotic
expansions near various values of the shape and argument parameters. We report
several numerical examples illustrating the efficiency of our methods. The
resulting code has been made available online.
