[page 639, §1]
valid for all , where . [639.1.3] Here denotes the largest integer smaller than . [639.1.4] Because of the recursion we consider only the case in the following.
[639.2.1] For the complex -plane is partitioned into four regions , , , . [639.2.2] In each region a different method is used for the calculation of the Mittag-Leffler function. [639.2.3] The central region
is the closure of the open disk of radius centered at the origin. [639.2.4] In our numerical calculations we have chosen throughout. [639.2.5] The regions are defined by
is the open wedge with opening angle . [639.2.6] For future reference we define also the (left) wedge
[639.2.7] Finally the region is defined as
where will be defined below.
[639.3.1] For , i.e. for the central disk region (with ), the Mittag-Leffler function is computed from truncating its Taylor series (1). [639.3.2] For given and accuracy we choose [page 640, §0] such that
holds. [640.0.1] The dependence of on the accuracy and other parameters is found as
[640.1.1] We start from the basic integral representation [2, p. 210]
where the path of integration in the complex plane starts and ends at and encircles the circular disc in the positive sense. [640.1.2] When this integral is evaluated several cases arise. [640.1.3] For we distinguish the cases and and compute the Mittag-Leffler function from the integral representations 
[640.1.4] For the integral representations read
[page 641, §1]
[641.1.1] The integrand is oscillatory but bounded over the integration interval. [641.1.2] Thus the integrals over can be evaluated numerically using any appropriate quadrature formula. [641.1.3] We use a robust Gauss-Lobatto quadrature.
[641.2.1] The integrals over involve infinite intervals. [641.2.2] For given (resp. ) and accuracy we approximate the integrals by truncation. [641.2.3] The error
depends on the truncation point . [641.2.4] For we find
while for we have
[641.3.1] The asymptotic expansions given in eqs. (21) and (22) on page 210 in Ref.  indicate that the Mittag-Leffler function exhibits a Stokes phenomenon. [641.3.2] The Stokes lines are the rays and it was recently shown in  that a Berry-type smoothing applies. [641.3.3] The exponentially improved asymptotic expansions will be used here for computing for . [641.3.4] More precisely for we use the following exponentially improved uniform asymptotic expansions  :
[page 642, §0] for with and , and
to truncate the asymptotic series. [642.0.5] Choosing for the machine precision we fix the radius in at