Code:
/* polylog.c
*
* Polylogarithms
*
*
*
* SYNOPSIS:
*/
double x, y, polylog();
int n;
y = polylog( n, x );
/* The polylogarithm of order n is defined by the series
*
* inf k
*  x
* Li (x) = >  .
* n  n
* k=1 k
*
*
* For x = 1,
*
* inf
*  1
* Li (1) = >  = Riemann zeta function (n) .
* n  n
* k=1 k
*
*
* When n = 2, the function is the dilogarithm, related to Spence's integral:
*
* x 1x
*  
*   ln(1t)   ln t
* Li (x) =   dt =   dt = spence(1x) .
* 2   t   1  t
*  
* 0 1
*
*
* See also the program cpolylog.c for the complex polylogarithm,
* whose definition is extended to x > 1.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain n # trials peak rms
* IEEE 0, 1 2 50000 6.2e16 8.0e17
* IEEE 0, 1 3 100000 2.5e16 6.6e17
* IEEE 0, 1 4 30000 1.7e16 4.9e17
* IEEE 0, 1 5 30000 5.1e16 7.8e17
*
*/
#include "cmath"
extern double PI;
/* polylog(4, 1x) = zeta(4)  x zeta(3) + x^2 A4(x)/B4(x)
0 <= x <= 0.125
Theoretical peak absolute error 4.5e18 */
#if UNK
static double A4[13] = {
3.056144922089490701751E2,
3.243086484162581557457E1,
2.877847281461875922565E1,
7.091267785886180663385E2,
6.466460072456621248630E3,
2.450233019296542883275E4,
4.031655364627704957049E6,
2.884169163909467997099E8,
8.680067002466594858347E11,
1.025983405866370985438E13,
4.233468313538272640380E17,
4.959422035066206902317E21,
1.059365867585275714599E25,
};
static double B4[12] = {
/* 1.000000000000000000000E0, */
2.821262403600310974875E0,
1.780221124881327022033E0,
3.778888211867875721773E1,
3.193887040074337940323E2,
1.161252418498096498304E3,
1.867362374829870620091E5,
1.319022779715294371091E7,
3.942755256555603046095E10,
4.644326968986396928092E13,
1.913336021014307074861E16,
2.240041814626069927477E20,
4.784036597230791011855E25,
};
#endif
#if DEC
static short A4[52] = {
0036772,0056001,0016601,0164507,
0037646,0005710,0076603,0176456,
0037623,0054205,0013532,0026476,
0037221,0035252,0101064,0065407,
0036323,0162231,0042033,0107244,
0035200,0073170,0106141,0136543,
0033607,0043647,0163672,0055340,
0031767,0137614,0173376,0072313,
0027676,0160156,0161276,0034203,
0025347,0003752,0123106,0064266,
0022503,0035770,0160173,0177501,
0017273,0056226,0033704,0132530,
0013403,0022244,0175205,0052161,
};
static short B4[48] = {
/*0040200,0000000,0000000,0000000, */
0040464,0107620,0027471,0071672,
0040343,0157111,0025601,0137255,
0037701,0075244,0140412,0160220,
0037002,0151125,0036572,0057163,
0035630,0032452,0050727,0161653,
0034234,0122515,0034323,0172615,
0032415,0120405,0123660,0003160,
0030330,0140530,0161045,0150177,
0026002,0134747,0014542,0002510,
0023134,0113666,0035730,0035732,
0017723,0110343,0041217,0007764,
0014024,0007412,0175575,0160230,
};
#endif
#if IBMPC
static short A4[52] = {
0x3d29,0x23b0,0x4b80,0x3f9f,
0x7fa6,0x0fb0,0xc179,0x3fd4,
0x45a8,0xa2eb,0x6b10,0x3fd2,
0x8d61,0x5046,0x2755,0x3fb2,
0x71d4,0x2883,0x7c93,0x3f7a,
0x37ac,0x118c,0x0ecf,0x3f30,
0x4b5c,0xfcf7,0xe8f4,0x3ed0,
0xce99,0x9edf,0xf7f1,0x3e5e,
0xc710,0xdc57,0xdc0d,0x3dd7,
0xcd17,0x54c8,0xe0fd,0x3d3c,
0x7fe8,0x1c0f,0x677f,0x3c88,
0x96ab,0xc6f8,0x6b92,0x3bb7,
0xaa8e,0x9f50,0x6494,0x3ac0,
};
static short B4[48] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x2e77,0x05e7,0x91f2,0x4006,
0x37d6,0x2570,0x7bc9,0x3ffc,
0x5c12,0x9821,0x2f54,0x3fd8,
0x4bce,0xa7af,0x5a4a,0x3fa0,
0xfc75,0x4a3a,0x06a5,0x3f53,
0x7eb2,0xa71a,0x94a9,0x3ef3,
0x00ce,0xb4f6,0xb420,0x3e81,
0xba10,0x1c44,0x182b,0x3dfb,
0x40a9,0xe32c,0x573c,0x3d60,
0x077b,0xc77b,0x92f6,0x3cab,
0xe1fe,0x6851,0x721c,0x3bda,
0xbc13,0x5f6f,0x81e1,0x3ae2,
};
#endif
#if MIEEE
static short A4[52] = {
0x3f9f,0x4b80,0x23b0,0x3d29,
0x3fd4,0xc179,0x0fb0,0x7fa6,
0x3fd2,0x6b10,0xa2eb,0x45a8,
0x3fb2,0x2755,0x5046,0x8d61,
0x3f7a,0x7c93,0x2883,0x71d4,
0x3f30,0x0ecf,0x118c,0x37ac,
0x3ed0,0xe8f4,0xfcf7,0x4b5c,
0x3e5e,0xf7f1,0x9edf,0xce99,
0x3dd7,0xdc0d,0xdc57,0xc710,
0x3d3c,0xe0fd,0x54c8,0xcd17,
0x3c88,0x677f,0x1c0f,0x7fe8,
0x3bb7,0x6b92,0xc6f8,0x96ab,
0x3ac0,0x6494,0x9f50,0xaa8e,
};
static short B4[48] = {
/*0x3ff0,0x0000,0x0000,0x0000,*/
0x4006,0x91f2,0x05e7,0x2e77,
0x3ffc,0x7bc9,0x2570,0x37d6,
0x3fd8,0x2f54,0x9821,0x5c12,
0x3fa0,0x5a4a,0xa7af,0x4bce,
0x3f53,0x06a5,0x4a3a,0xfc75,
0x3ef3,0x94a9,0xa71a,0x7eb2,
0x3e81,0xb420,0xb4f6,0x00ce,
0x3dfb,0x182b,0x1c44,0xba10,
0x3d60,0x573c,0xe32c,0x40a9,
0x3cab,0x92f6,0xc77b,0x077b,
0x3bda,0x721c,0x6851,0xe1fe,
0x3ae2,0x81e1,0x5f6f,0xbc13,
};
#endif
#ifdef ANSIPROT
extern double spence ( double );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double zetac ( double );
extern double md_pow ( double, double );
extern double md_powi ( double, int );
extern double md_log ( double );
extern double fac ( int i );
extern double md_fabs (double);
double polylog (int, double);
#else
extern double spence(), polevl(), p1evl(), zetac();
extern double md_pow(), md_powi(), md_log();
extern double fac(); /* factorial */
extern double md_fabs();
double polylog();
#endif
extern double MACHEP;
double
polylog (n, x)
int n;
double x;
{
double h, k, p, s, t, u, xc, z;
int i, j;
/* This recurrence provides formulas for n < 2.
d 1
 Li (x) =  Li (x) .
dx n x n1
*/
if (n == 1)
{
p = 1.0  x;
u = x / p;
s = u * u + u;
return s;
}
if (n == 0)
{
s = x / (1.0  x);
return s;
}
/* Not implemented for n < 1.
Not defined for x > 1. Use cpolylog if you need that. */
if (x > 1.0  n < 1)
{
mtherr("polylog", DOMAIN);
return 0.0;
}
if (n == 1)
{
s = md_log (1.0  x);
return s;
}
/* Argument +1 */
if (x == 1.0 && n > 1)
{
s = zetac ((double) n) + 1.0;
return s;
}
/* Argument 1.
1n
Li (z) =  (1  2 ) Li (z)
n n
*/
if (x == 1.0 && n > 1)
{
/* Li_n(1) = zeta(n) */
s = zetac ((double) n) + 1.0;
s = s * (md_powi (2.0, 1  n)  1.0);
return s;
}
/* Inversion formula:
* [n/2] n2r
* n 1 n  md_log (z)
* Li (z) + (1) Li (1/z) =   md_log (z) + 2 >  Li (1)
* n n n!  (n  2r)! 2r
* r=1
*/
if (x < 1.0 && n > 1)
{
double q, w;
int r;
w = md_log (x);
s = 0.0;
for (r = 1; r <= n / 2; r++)
{
j = 2 * r;
p = polylog (j, 1.0);
j = n  j;
if (j == 0)
{
s = s + p;
break;
}
q = (double) j;
q = md_pow (w, q) * p / fac (j);
s = s + q;
}
s = 2.0 * s;
q = polylog (n, 1.0 / x);
if (n & 1)
q = q;
s = s  q;
s = s  md_pow (w, (double) n) / fac (n);
return s;
}
if (n == 2)
{
if (x < 0.0  x > 1.0)
return (spence (1.0  x));
}
/* The power series converges slowly when x is near 1. For n = 3, this
identity helps:
Li (x/(1x)) + Li (1x) + Li (x)
3 3 3
2 2 3
= Li (1) + (pi /6) md_log(1x)  (1/2) md_log(x) md_log (1x) + (1/6) md_log (1x)
3
*/
if (n == 3)
{
p = x * x * x;
if (x > 0.8)
{
u = md_log(x);
s = u * u * u / 6.0;
xc = 1.0  x;
s = s  0.5 * u * u * md_log(xc);
s = s + PI * PI * u / 6.0;
s = s  polylog (3, xc/x);
s = s  polylog (3, xc);
s = s + zetac(3.0);
s = s + 1.0;
return s;
}
/* Power series */
t = p / 27.0;
t = t + .125 * x * x;
t = t + x;
s = 0.0;
k = 4.0;
do
{
p = p * x;
h = p / (k * k * k);
s = s + h;
k += 1.0;
}
while (md_fabs(h/s) > 1.1e16);
return (s + t);
}
if (n == 4)
{
if (x >= 0.875)
{
u = 1.0  x;
s = polevl(u, A4, 12) / p1evl(u, B4, 12);
s = s * u * u  1.202056903159594285400 * u;
s += 1.0823232337111381915160;
return s;
}
goto pseries;
}
if (x < 0.75)
goto pseries;
/* This expansion in powers of md_log(x) is especially useful when
x is near 1.
See also the pari gp calculator.
inf j
 z(nj) (md_log(x))
polylog(n,x) = > 
 j!
j=0
where
z(j) = Riemann zeta function (j), j != 1
n1

z(1) = md_log(md_log(x)) + > 1/k

k=1
*/
z = md_log(x);
h = md_log(z);
for (i = 1; i < n; i++)
h = h + 1.0/i;
p = 1.0;
s = zetac((double)n) + 1.0;
for (j=1; j<=n+1; j++)
{
p = p * z / j;
if (j == n1)
s = s + h * p;
else
s = s + (zetac((double)(nj)) + 1.0) * p;
}
j = n + 3;
z = z * z;
for(;;)
{
p = p * z / ((j1)*j);
h = (zetac((double)(nj)) + 1.0);
h = h * p;
s = s + h;
if (md_fabs(h/s) < MACHEP)
break;
j += 2;
}
return s;
pseries:
p = x * x * x;
k = 3.0;
s = 0.0;
do
{
p = p * x;
k += 1.0;
h = p / md_powi(k, n);
s = s + h;
}
while (md_fabs(h/s) > MACHEP);
s += x * x * x / md_powi(3.0,n);
s += x * x / md_powi(2.0,n);
s += x;
return s;
}