# How to use polylogarithm function in c++?

Is there any preprocessor directives that could be used to use the polylog function? Or is it included in cmath? If so, do you call it by Li or by polylog?

EDIT: What I really am trying to do is give an analytical value for the indefinite integral of the function

$$\frac{x^3}{e^{x} -1}$$

which involves polylogarithm functions. But if anyone has a suggestion for another way to integrate this function analytically I'd be welcome to any ideas.

• For searching purposes: the function being considered in the OP is related to the Debye functions. This note might be of some of use. May 8, 2013 at 12:57
• Another close relation is the incomplete Fermi-Dirac integral. May 15, 2013 at 16:35

There's a GPL'd C library, ANANT - Algorithms in Analytic Number Theory by Linas Vepstas, which includes multiprecision implementation of the polylogarithm, building on GMP.

From its README file:

This project contains ad-hoc implementations of assorted analytic functions of interest in number theory, including the gamma function, the Riemann zeta function, the polylogarithm, and the Minkowski question mark function. The implementation uses the Gnu Multi-Precision library (GMP) to perform all low-level operations. The code herein is licensed under the terms of the Gnu GPLv3 license.

The GSL (GNU Scientific Library) apparently only has the dilogarithm function. However following a hint from @J.M. one finds the Debye function which gives the ulterior integral (up to a scalar multiple) implemented in double precision (see GSL 7.10 Debye Functions orders 1 through 6):

$$D_n(x) = \frac{n}{x^n} \int_0^x \frac{t^n dt}{e^t - 1}$$

Symbolic integration software such as Mathematica or Maxima provides:

$$\int_0^x \frac{t^3 dt}{e^t-1} = 6 \operatorname{Li}_4(e^x) - 6x \operatorname{Li}_3(e^x) + 3x^2 \operatorname{Li}_2(e^x) + x^3 \log(1-e^x) - \frac{x^4}{4} - \frac{\pi^4}{15}$$

The left-hand side is obviously a purely real value if $$x \gt 0$$, but the polylogarithms shown will be complex-valued (because $$e^x > 1$$, and so equality hinges on total cancellation of imaginary parts). We can avoid the need for complex arithmetic in this case by substituting the expression:

$$\int_0^x \frac{t^3 dt}{e^t-1} = -6 \operatorname{Li}_4(e^{-x}) - 6x \operatorname{Li}_3(e^{-x}) - 3x^2 \operatorname{Li}_2(e^{-x}) - x^3\operatorname{Li}_1(e^{-x}) + \frac{\pi^4}{15}$$

This is an improvement because with polylogarithm arguments in $$[0,1]$$, the results are purely real values. Note the proper result when $$x = 0$$ is zero, and this is achieved by cancellation between the leading term and the constant. Thus relative error could be an issue for small positive values of $$x$$.

Note that our mysterious constant $$\pi^4/15$$ is the limiting upper bound on these (monotone increasing) integrals:

$$\int_0^\infty \frac{t^3 dt}{e^t - 1} = \Gamma(4) \zeta(4) = 6 \cdot \frac{\pi^4}{90}$$

We can now revisit the title question, How to use polylogarithm function in c++? The point is worth making that there is no standard implementation of polylogarithm functions for C or even C++. If the goal is to avoid any additional library for your implementation, it pretty well sets you out to rolling your own routines, perhaps along the lines suggested by the David C. Wood paper that GertVdE’s Answer links to.

Besides the multiprecision routines suggested in the first part of my Answer, there is a mature (free) double precision math library in Cephes by Stephen L. Moshier which implements both real (polylog) and complex (cpolylog) versions of the polylogarithm special functions. Although their accuracy depends in part on the underlying standard mathematical functions of C, the Cephes source documentation reports tests and theoretical peak errors for orders 1 through 4 at about the limits of double precision.

Alternatively you may wish to use other software to directly check (not referencing polylogarithms) the quadrature routines you wrote for your integral. As I sketch out in this Math.SE Question, the power series centered at the origin for the integral has limited convergence, but this can be mitigated by using a continued fraction expansion instead.

For immediate gratification I recommend the (free) numerical quadrature QUADPACK routines included in Maxima, specifically quad_qag. For example find the integral over [0,5] with this Maxima command:

(%i1) quad_qag(x^3/(%e^x - 1), x, 0, 5, 2);
(%o1) [4.899892158330582,5.4399730923588665*10^-14,21,0]


Of the input arguments only the last one bears an explanation. The fifth argument to quad_qag specifies what rule to apply in adaptive quadrature. Possible values are 1 to 6, and give increasing sophistication/accuracy. The output line gives first the numerical quadrature, followed by an estimate of its absolute error, the number of subintervals/steps used, and a return code (here zero means no error or special conditions found).

First of all, you should choose based on your application if you need high precision arithmetic (i.e. will you be happy with just IEEE double precision results for the polylog functions or do you need higher precision)? If you do need high precision, you can look in the family of tools around the GMP library.

If you don't, you can use approximations. Some literature research pointed me to this article. At the end of the article, is a "selection table": based on the arguments of the polylogs that you need, you can select an approximation formula. But be careful to check stability and accuracy.

If you don't need too many evaluations (not in a nested loop), I would just go for numerical quadrature using the double exponentional method.

• Precision needed would be around 10^-6. Do you know any other way which doesn't involve using any extra libraries? Apr 29, 2013 at 10:05
• @flamingohats: then I would go for the expressions in the paper I linked to in my answer. What is your range for $x$? Apr 29, 2013 at 15:08
• @flamingohats: Are you saying you need to approximate $\int_0^x \frac{t^3}{e^t - 1} dt$ for $x \in [0,10]$ or something of that nature? I don't follow why "there will only be two evaluations needed", unless you mean that this is an additional requirement/ideal for the quadrature or other approximation. I'm thinking we can determine a polynomial or rational approximation that gives your required accuracy. Apr 30, 2013 at 15:07
• sorry I got confused. I have approximated the integral from 0.65 to 5.025 using the Trapezoidal method, and I need a formula to find the exact value so that I can compare the approximation with an analytical value. I know this will be approximated because it is a floating point number, so a precision of 1e-6 will be fine. If I can somehow learn how to input the polylog function into the IDE it should work. May 1, 2013 at 16:36
• @flamingohats: it still is not clear whether you need a one shot value ($\int_{0.65}^{5.205} \frac{t^3}{e^t-1}dt = 4.8498308528256668370925$) or if you really need it programmed in your code. If it is the first, why not use Sage, Euler Toolbox, ... to check a number of crucial values for you (either using their polylog implementation or their sophisticated quadrature rules)? If it is the second, implement a "traditional" quadrature method (like you did) and again use Sage, Euler, ... to test and validate it. May 2, 2013 at 5:55

The GSL also has complete Fermi-Dirac integrals, $F_j(x)$ for integer $j$ and $j=-\frac12, \frac12, \frac32$. These functions are equivalent to polylogs $$F_j(x) = -\text {Li}_{j+1}(-e^{-x})$$ Though note the restriction to negative polylog arguments for real $x$.