A fast and simple algorithm for the computation of the Lerch transcendent

TL;DR

This paper presents a fast, reliable algorithm using Gauss-Laguerre quadrature for computing the Lerch transcendent, including error estimates and a MATLAB implementation, enabling arbitrary precision calculations.

Contribution

It introduces a novel, efficient algorithm for Lerch transcendent computation with error control and practical MATLAB code, improving upon existing methods.

Findings

Error estimates enable precise quadrature node selection

Algorithm achieves high accuracy with fewer nodes

Numerical tests confirm reliability and efficiency

Abstract

This paper deals with the computation of the Lerch transcendent by means of the Gauss-Laguerre formula. An a priori estimate of the quadrature error, that allows to compute the number of quadrature nodes necessary to achieve an arbitrary precision, is derived. Exploiting the properties of the Gauss-Laguerre rule and the error estimate, a truncated approach is also considered. The algorithm used and its Matlab implementation are reported. The numerical examples confirm the reliability of this approach.

Tables

Table 1. Table 1: Results of Algorithm 1 for the polylogarithm Li s ​ ( r ​ e i ​ τ ​ π ) subscript Li 𝑠 𝑟 superscript 𝑒 𝑖 𝜏 𝜋 \text{Li}_{s}\left(re^{i\tau\pi}\right) .

Table 2. Table 2: Results of Algorithm 1 for the Dirichlet beta function β ​ ( s ) 𝛽 𝑠 \beta(s) .

Table 3. Table 3: Results of Algorithm 1 for the Dirichlet eta function η ​ ( s ) 𝜂 𝑠 \eta(s) .

Table 4. Table 4: Results of Algorithm 1 for general cases of the Lerch transcendent Φ ​ ( r ​ e i ​ τ ​ π , s , a ) Φ 𝑟 superscript 𝑒 𝑖 𝜏 𝜋 𝑠 𝑎 \Phi\left(re^{i\tau\pi},s,a\right) .

Equations

$$\Phi (z,s,a)=\sum _{j=0}^{\infty }\frac{z^j}{(j+a{)}^s},$$

$$\Phi (z,s,a)=\frac{1}{\Gamma (s)}{\int }_0^{\infty }\frac{x^{s-1}{e}^{-ax}}{1-z{e}^{-x}}dx,z\in \mathbb{C}\setminus [1,+\infty ),\Re (s)>0,\Re (a)>0,$$

$$\Phi (z,s,a)=\frac{1}{\Gamma (s)a^s}I(z),$$

$$I(z)={\int }_0^{\infty }\frac{t^{s-1}{e}^{-t}}{1-z{e}^{-\frac{t}{a}}}dt.$$

$$I_n(z)=\sum _{j=1}^n{w}_j^{(n)}{f}_z\left(t_j^{(n)}\right),f_z(t)=\frac{1}{1-z{e}^{-\frac{t}{a}}},$$

$$\Phi (z,s,a)\approx {\Phi }_n(z,s,a)=\frac{1}{\Gamma (s)a^s}{I}_n(z).$$

$$E_n=\Phi (z,s,a)-{\Phi }_n(z,s,a)=\frac{1}{\Gamma (s)a^s}{e}_n,$$

$$e_n=I(z)-I_n(z),$$

$$t_k=a(\ln |z|+i(\arg (z)+2k\pi )),-\pi \le \arg (z)\le \pi ,k\in \mathbb{Z}.$$

$$e_n=\frac{1}{2\pi i}{\int }_{\Gamma }\frac{q_n(w)}{L_n^{(s-1)}(w)}{f}_z(w)dw,$$

$$\Re \sqrt{-w}=\ln R.$$

$$w=\frac{y^2}{4{\left(\ln R\right)}^2}-{\left(\ln R\right)}^2+iy,y\in \mathbb{R}.$$

$$\Gamma ={\Gamma }_R\bigcup \left(\bigcup _{k=-N}^N{C}_k\right),$$

$$\frac{q_n(w)}{L_n^{(s-1)}(w)}\sim -2e^{-i\pi (s-1)}{w}^{s-1}{e}^{-w}\frac{K_{s-1}\left(2\sqrt{m}{e}^{-i\pi /2}\sqrt{w}\right)}{I_{s-1}\left(2\sqrt{m}{e}^{-i\pi /2}\sqrt{w}\right)},$$

$$I_{\alpha }(t)$$

$$K_{\alpha }(t)$$

$$\frac{K_{\alpha }(t)}{I_{\alpha }(t)}\sim \pi e^{-2t}\left(1+\mathcal{O}\left(\frac{1}{t}\right)\right).$$

$$\frac{q_n(w)}{L_n^{(s-1)}(w)}\sim 2\pi e^{-is\pi }{w}^{s-1}{e}^{-w}{e}^{-4\sqrt{m}\sqrt{-w}},n\to +\infty ,$$

$$e_n\sim -i{e}^{-is\pi }{\int }_{\Gamma }{w}^{s-1}{e}^{-w}{e}^{-4\sqrt{m}\sqrt{-w}}{f}_z(w)dw.$$

$$g(w):=w^{s-1}{e}^{-w}{e}^{-4\sqrt{m}\sqrt{-w}},$$

$$e_n\sim -i{e}^{-is\pi }\left[{\int }_{{\Gamma }_R}g(w)f_z(w)dw-2\pi i\sum _{k=-N}^N\mathrm{Res}\left(g(w)f_z(w),t_k\right)\right],$$

$${\int }_{{\Gamma }_R}g(w)f_z(w)dw$$

$$\chi (w)=e^{-4i\sqrt{m}\Im \sqrt{-w}}.$$

$$\mathrm{Res}\left(g(w)f_z(w),t_k\right)$$

$${\ln }_k(z)=\ln |z|+i(\arg (z)+2k\pi ),k\in \mathbb{Z},$$

$$\begin{array}{rl}{e}_n& \sim -i{e}^{-is\pi }[R^{-4\sqrt{m}}{\int }_{{\Gamma }_R}{w}^{s-1}{e}^{-w}\chi (w)f_z(w)dw\\ & -2\pi i{a}^s{z}^{-a}\sum _{k=-N}^N{\left[{\ln }_k(z)\right]}^{s-1}\chi (t_k)R_k^{-4\sqrt{m}}].\end{array}$$

$$e_n\sim -2\pi e^{-is\pi }{a}^s{\left[{\ln }_0(z)\right]}^{s-1}{z}^{-a}\chi (t_0)R_0^{-4\sqrt{m}}.$$

$$e_n$$

$$|e_n|$$

$$=4\pi a^s|z{|}^{-a}|{\ln }_0(z){|}^{s-1}{e}^{-4\sqrt{m}\Re \left(\sqrt{-t_0}\right)}=:{ϵ}_n,$$

$$|E_n|\approx \frac{4\pi }{\Gamma (s)}|z{|}^{-a}|{\ln }_0(z){|}^{s-1}{e}^{-4\sqrt{m}\Re \left(\sqrt{-t_0}\right)}=:{\mathcal{E}}_n.$$

Full text

1 Introduction

In this work we consider the Lerch transcendent introduced in [10] and defined as (see [14, p.628 n.25.14.1])

[TABLE]

where $a\neq 0,-1,-2,\ldots$, $|z|<1$ or $\Re(s)>1$, $|z|=1$. For other values of $z$, $\Phi(z,s,a)$ is defined by analytic continuation. For a recent investigation on the analytic properties of $\Phi(z,s,a)$, we refer to [11].

The Lerch transcendent and its special cases, such as the polylogarithm (corresponding to $a=1$), appears for instance in quantum Bose and Fermi statistics. In particular, the particle number density, the pressure, the internal energy and the entropy of the ideal quantum gases of Fermi and Bose can be expressed in terms of the polylogarithm [5]. Another application of the Lerch transcendent is in biophysics. Indeed, some discrete distributions, related to the Lerch transcendent by the normalization constant of the associated probability mass functions, are used to establish the statistical composition of DNA and protein sequences [2].

In order to evaluate the Lerch transcendent, some authors have used series representations and asymptotic expansions. Among the others, we quote [6, 3] for an overview of properties, identities and numerical methods for the computation of the Lerch transcendent. Moreover, in [2], combined non-linear condensation transformation, which is an algorithm that allows to evaluate slowly convergent nonalternating series, is considered. More recently, in [13], starting from the Hermite-type integral representation of the Lerch transcendent, a new uniform asymptotic expansion of $\Phi(z,s,a)$, for large order of the parameters $a,s$ and argument $z$, is derived.

In this work, exploiting the integral representation

[TABLE]

where $\Gamma$ is the Gamma function, we compute $\Phi(z,s,a)$ by employing the Gauss-Laguerre formula and, starting from the analysis given in [4], we derive an accurate error estimate for the quadrature error. To the best of our knowledge, the Gauss-Laguerre formula has never been used to this purpose, even if the generalized Laguerre weight clearly appears in the integral representation (1). Moreover, after removing the weight, what remains is bounded, and therefore the use of the Gauss-Laguerre method appears a natural choice for this problem. Finally, since the weights of this rule decay exponentially, we present a reliable algorithm which allows to reduce the number of function evaluations, without significant loss of accuracy in the computation of the integral.

Throughout this work we use the symbol $\approx$ to indicate a generic approximation. Whereas, the symbol $\sim$ is used to express the asymptotic equality.

The paper is organized as follows. In Section 2, working with $s$ and $a$ real, we present the Gauss-Laguerre approach, together with the error analysis and some numerical experiments. In Section 3 we introduce the truncated approach and the algorithm that allows to compute the Lerch transcendent with a prescribed accuracy. Some tables with experiments are given in Section 4. In Section 5 we show how to extend the method for $s$ and $a$ complex. An appendix contains an auxiliary result and the Matlab code we have used.

2 The Gauss-Laguerre approach

For simplicity we assume $s,a\in\mathbb{R}$. In order to employ the Gauss-Laguerre rule, we consider the change of variable $ax=t$ in (1), that leads to

[TABLE]

where

[TABLE]

Denoting by

[TABLE]

the corresponding $n$-point Gauss-Laguerre rule, where $t_{j}^{(n)}$ and $w_{j}^{(n)}$, $j=1,\ldots,n$, are the nodes and the weights respectively, we obtain the approximation

[TABLE]

2.1 Error analysis

In order to derive an estimate for the error

[TABLE]

with

[TABLE]

we first need to locate the poles of the function $f_{z}$ defined in (2). By solving $1-ze^{-t/a}=0$, we obtain

[TABLE]

It is known that for the error of the Gauss-Laguerre formula it holds

[TABLE]

where $L_{n}^{(s-1)}$ is $n$-th degree generalized Laguerre polynomial and $q_{n}$ is the corresponding associated function (see [7, 15]). $\Gamma$ is a contour in the complex plane containing the set $[0,+\infty)$ but no singularity of the function $f_{z}$ can lie on or within $\Gamma$. For a suitable choice of $\Gamma$, we consider, for $R>1$, the parabola of the complex plane defined by the equation

[TABLE]

It is immediate to verify that the parabola, which we denote by $\Gamma_{R}$, can be rewritten as

[TABLE]

It is symmetric with respect to the real axis, with vertex in $-\left(\ln R\right)^{2}$ and convexity oriented towards the positive real axis. The parabola degenerates to $[0,+\infty)$ as $R\rightarrow 1$.

Let $C_{k}$ be an arbitrary small circle surrounding the pole $t_{k}$. Assuming $-\pi<\arg(z)<\pi$, the pole closest to the real axis is $t_{0}$ (see (5)), and therefore for any non negative integer $N$ we can define

[TABLE]

with $R$ such that the parabola $\Gamma_{R}$ contains in its interior the poles $t_{-N},...,t_{N}$ but none of the others.

In order to evaluate (6), with $\Gamma$ as in (9), we first observe that, for n $\rightarrow+\infty$, [7, p.33]

[TABLE]

where $m=n+\frac{s}{2}$ and $I$, $K$ are the modified Bessel functions of the first and second kind, respectively. Since $0<\arg(w)<2\pi$, we can use for $\alpha>-1$ the asymptotic expansions (see [1, p.377 n.9.7.1-2])

[TABLE]

valid for large $|t|$, $\left|\arg(t)\right|<\frac{\pi}{2}$, and therefore,

[TABLE]

By using (11) in (10), we find

[TABLE]

and therefore

[TABLE]

By defining for simplicity

[TABLE]

and with $\Gamma$ as in (9), we obtain

[TABLE]

where ${\rm Res(\cdot,\cdot)}$ denotes the residue. At this point, by (12) and (7)

[TABLE]

where

[TABLE]

We remark that, since $|f_{z}|$ is bounded, the above integral is also bounded. Moreover, since ${\rm Res}\left(f_{z}(w),t_{k}\right)=a$, $k\in\mathbb{Z}$, by (5) we have

[TABLE]

where we use the notation

[TABLE]

and $R_{k}<R$ (for $-N\leq k\leq N$) represents the parabola passing through $t_{k}$, that is, $R_{k}$ solves $\Re\sqrt{-t_{k}}=\ln R_{k}$. Joining the above results we finally obtain

[TABLE]

Under the assumption $\left|\arg(z)\right|<\pi$, we have $R_{0}<R_{k}<R$, for $|k|=1,...,N$, and therefore, by collecting the term involving $R_{0}$,

[TABLE]

This situation in shown in Figure 1a, where $z=\sqrt{2}e^{i\frac{\pi}{4}}$.

In the case of $\arg(z)=\pm\pi$ we have that the poles are symmetric with respect to the real axes and it holds $\overline{t_{k}}=t_{-k-1}$ if $\arg(z)=\pi$, whereas $\overline{t_{-k}}=t_{k+1}$ if $\arg(z)=-\pi$, $k\geq 0$ (see Figure 1b, in which $z=e^{i\pi}$). In both situations, formula (15) needs to be replaced by the sum of the terms corresponding to the poles $t_{0}$ and $t_{-1}$, for $\arg(z)=\pi$, and $t_{0}$ and $t_{1}$ for $\arg(z)=-\pi$. In both cases, we obtain the same result. Assuming for instance $\arg(z)=\pi$, we have that $\overline{t_{0}}=t_{-1}$, $R_{0}=R_{-1}$ and therefore

[TABLE]

Working with the modulus, the only essential difference between the above formula and (15) is a factor 2. Moreover, since the latter should also be used for $\arg(z)$ close to $\pm\pi$, independently of $z$ we finally consider the error estimate

[TABLE]

for the integral, and therefore (see (3))

[TABLE]

2.2 Some numerical experiments

In Figure 2 we show the quality of estimate (17) on some examples. Specifically, we consider the polylogarithm [14, n.25.12] (Figure 2a)

[TABLE]

the Dirichlet beta function [1, p.807, n.23.2.21] (Figure 2b)

[TABLE]

and the Dirichlet eta function [1, p.807, n.23.2.19] (Figure 2c)

[TABLE]

These functions can be written in terms of $\Phi(z,s,a)$ for particular values of the parameters, that is,

[TABLE]

We remark that the value $s=3/2$ chosen for the polylogarithm represents the case of the three-dimensional Bose gas in a box. Finally, in Figure 2d we consider a general situation. As for the computation of the nodes and weights of the generalized Gauss-Laguerre rule (c.f. (2)), we employ the Matlab routine $\mathtt{lagpts.m}$ [16], based on the traditional Golub-Welsch algorithm (see [8, 7]). We point out that the oscillations of the error in Figure 2a-b-c are due to the term $\chi(t_{0})$ (cf. (14)- (15)).

3 A truncated approach

Having at disposal an accurate error estimate, in this section we derive a reliable algorithm that allows to reduce the number of evaluations of the function $f_{z}$ (see (2)), required by the quadrature rule. This is possible because the weights of the Gauss-Laguerre rule decay exponentially and, moreover, because the function $f_{z}$ is bounded. Indeed, a simple analysis, reported in Appendix, shows that

[TABLE]

Now, let $\tau_{n}$ be the solution of

[TABLE]

where $\epsilon_{n}$ is defined in (2.1). By using the relation [9, p.942 n.8.357]

[TABLE]

we approximate $\tau_{n}$ by solving with respect to $x$ the equation

[TABLE]

The solution is given by

[TABLE]

where $W$ denotes the Lambert W function. In particular, we obtain

[TABLE]

for $0<s<1$, and

[TABLE]

for $s>1$. It is known that ([14, n. 4.13.10-11])

[TABLE]

and therefore $x\sim g_{n}(s)$, where

[TABLE]

Note that $g_{n}(s)$ has a removable singularity in $s=1$, so that we can define

[TABLE]

that is also the exact solution of (19), with respect to $\tau_{n}$, for $s=1$. In conclusion, we have that $\tau_{n}\sim g_{n}(s)$, for $s>0$.

We are now on the point to introduce the truncated rule

[TABLE]

where $k_{n}$ is the smallest integer such that $x_{j}^{(n)}\geq g_{n}(s)$, for each $j\geq k_{n}$. As for the error, we have

[TABLE]

Now, using the bound [12, eqs. 2.4-2.7]

[TABLE]

where $c$ is a constant independent of $j,n$ and close to $1$ for large $n$, we have

[TABLE]

Finally, for the truncated rule we obtain the estimate

[TABLE]

In order to understand the decay rate with respect to $k_{n}$ we use the relations [1, n. 22.16.8]

[TABLE]

and [1, n. 9.5.12]

[TABLE]

to finally obtain

[TABLE]

At this point, it is possible to derive an analytical approximation of $k_{n}$ by imposing $x_{k}^{(n)}=g_{n}(s)$ with the help of (25). In order to detect the new asymptotic behavior, we further simplify the computation by neglecting the term involving $s$ in (22), that is, we solve with respect to $k$

[TABLE]

We remark, however, that this simplification is not used in the algorithm presented below. By (2.1) and defining

[TABLE]

we find

[TABLE]

so that, for $n$ large enough,

[TABLE]

By using this relation in (2.1) we finally obtain

[TABLE]

which expresses the speed up attainable with the truncation. Indeed, the number of function evaluations is now raised to the power of $2/3$, whereas in (2.1) the power is $1/2$. This decay rate is well visible in Figure 3, where we show the benefits of the truncated rule on some examples. In particular, as in Section 2.2, we consider the polylogarithm (Figure 3a), the Dirichlet beta function (Figure 3b), the Dirichlet eta function (Figure 3c) and a general situation (Figure 3d).

We summarize the basic steps for the computation of the Lerch transcendent by using the truncated Laguerre rule and with a prescribed error tolerance $\mathcal{E}$ in the following algorithm.

Algorithm 1** (Truncated Laguerre rule)**

Given $z,s,a,\mathcal{E}$

compute the corresponding tolerance for the integral $\epsilon=a^{s}\Gamma(s)\frac{\mathcal{E}}{2}$ (see (3) and (24)) 2. 2.

compute $K_{z}$ as in (18) 3. 3.

compute $R_{0}$ such that $t_{0}\in\Gamma_{R_{0}}$ by using (8) 4. 4.

solve $\epsilon_{n}=\epsilon$ with respect to $m$ (see (2.1)) and then set $n=m-s/2$ 5. 5.

estimate $k_{n}$ by solving, with respect to $k$, $x_{k}^{(n)}=g_{n}(s)$, by using (25) and (22) (in the implementation it is convenient to add one or two more points in order to prevent the under estimation of the error, because of the large number of approximations used) 6. 6.

compute the first $k_{n}$ nodes and weights $x_{j}^{(n)},w_{j}^{(n)}$, $j=1,...,k_{n}$, of the $n$-point Laguerre rule 7. 7.

approximate $\Phi(z,s,a)$ with $\Phi_{k_{n}}(z,s,a)=\frac{1}{\Gamma(s)a^{s}}I_{k_{n}}(z)$ (see (23))

4 Numerical experiments

In this section, working with prescribed error tolerances $\mathcal{E}=1e-10$ and $\mathcal{E}=1e-14$, we test Algorithm 1 on several examples. In particular, we present some tables in which, for each set of parameters, we report the corresponding values of $n$ and $k_{n}$, given by steps $4$-$5$, together with the final error obtained by using a reference solution. Specifically, we consider the polylogarithm $\text{Li}_{s}\left(re^{i\tau\pi}\right)$ in Table 1, the Dirichlet beta function $\beta(s)$ in Table 2, the Dirichlet eta function $\eta(s)$ in Table 3 and general cases of the Lerch transcendent $\Phi\left(re^{i\tau\pi},s,a\right)$ in Table 4, for different values of the parameters $r,\tau,s,a$. In Tables 1 and 4, $\tau$ is set in order to consider the arguments $\pi,\frac{3}{4}\pi,\frac{\pi}{2},\frac{\pi}{4}$. We can see that, except in few rare cases, the prescribed tolerance is achieved.

5 Complex case

Starting from the integral representation (1)

[TABLE]

by using the change of variable $x=\frac{t}{\Re(a)}$, we obtain

[TABLE]

The function

[TABLE]

inside the integral is the main difference with respect to the case of $s,a\in\mathbb{R}$. Some experiments have revealed that the case of $\Im(a)\neq 0$ does not constitute a problem for the Laguerre rule, unless $\Im(a)\gg\Re(a)$. On the other side, the case of $\Im(s)\neq 0$ may be difficult to handle. Indeed, as $t\rightarrow 0$, the real and the imaginary part of the function $\varphi(t)$ oscillate with increasing frequency and the Laguerre rule appears to be inadequate (Figure 4a). The situation is overtaken whenever $\Re(s)\gg\left|\Im(s)\right|$ and $\Re(s)>1$, since the term $t^{\Re(s)-1}$ makes the oscillations negligible near zero (see Figure 4b).

6 Conclusion

In this work we have employed the generalized Gauss-Laguerre formula to compute the Lerch transcendent $\Phi(z,s,a)$ for $z\in\mathbb{C}\setminus[1,+\infty)$ and $s,a>0$. We have derived sharp error estimates that enable to know a priori the number of quadrature points necessary to achieve a prescribed accuracy and to truncate the rule. We have tested the arising algorithm on several examples and the results confirm the reliability of this approach. The extension to $s$ and $a$ complex is not theoretically analyzed. Anyway, at a first glance the Laguerre rule appears robust if the imaginary parts are relatively small with respect to the moduli. On the contrary, the high frequency of the oscillations makes the method unsuited.

Appendix A Bounds for the integrand function

In order to derive the bound (18), we define $x=e^{-t/a}$, $t\in[0,+\infty)$, so that

[TABLE]

The function $g(x)=1+|z|^{2}x^{2}-2\Re(z)x$ is a parabola and its minimum in $[0,1]$ defines $K_{z}$. The vertex is at $x_{V}=\frac{\Re(z)}{|z|^{2}}$ and we have three cases. Let $D=\mathbb{C}\setminus[0,+\infty)$ be the domain of definition of integral (1) with respect to $z$. Let $A=\left\{z\in\mathbb{C}\mid\Re(z)\leq 0\right\}$, $B=\left\{z\in\mathbb{C}\mid\left|z-\frac{1}{2}\right|<\frac{1}{2}\right\}$ and $C=D\setminus(A\cup B)$. Note that $A,B,C$ are mutually disjoint and $A\cup B\cup C=D$. Moreover, the set $B$ is the solution of the inequality $\Re(z)>|z|^{2}$. At this point, if $z\in A$, then $x_{V}\leq 0$ and $g(x)\geq g(0)$, and therefore

[TABLE]

If $z\in B$, then $x_{V}>1$, $g(x)\geq g(1)$ and we find

[TABLE]

Finally, if $z\in C$, then $0<x_{V}\leq 1$, $g(x)\geq g(x_{V})$ and we obtain

[TABLE]

Appendix B Matlab code

Acknowledgements

This work was partially supported by GNCS-INdAM, FRA-University of Trieste and CINECA under HPC-TRES program award number 2019-04. The authors are member of the INdAM research group GNCS.