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.
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.
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 $.
