lib/libm/polylog.c

   1 /*                                                      polylog.c
   2  *
   3  *      Polylogarithms
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, y, polylog();
  10  * int n;
  11  *
  12  * y = polylog( n, x );
  13  *
  14  *
  15  * The polylogarithm of order n is defined by the series
  16  *
  17  *
  18  *              inf   k
  19  *               -   x
  20  *  Li (x)  =    >   ---  .
  21  *    n          -     n
  22  *              k=1   k
  23  *
  24  *
  25  *  For x = 1,
  26  *
  27  *               inf
  28  *                -    1
  29  *   Li (1)  =    >   ---   =  Riemann zeta function (n)  .
  30  *     n          -     n
  31  *               k=1   k
  32  *
  33  *
  34  *  When n = 2, the function is the dilogarithm, related to Spence's integral:
  35  *
  36  *                 x                      1-x
  37  *                 -                        -
  38  *                | |  -ln(1-t)            | |  ln t
  39  *   Li (x)  =    |    -------- dt    =    |    ------ dt    =   spence(1-x) .
  40  *     2        | |       t              | |    1 - t
  41  *               -                        -
  42  *                0                        1
  43  *
  44  *
  45  *  See also the program cpolylog.c for the complex polylogarithm,
  46  *  whose definition is extended to x > 1.
  47  *
  48  *  References:
  49  *
  50  *  Lewin, L., _Polylogarithms and Associated Functions_,
  51  *  North Holland, 1981.
  52  *
  53  *  Lewin, L., ed., _Structural Properties of Polylogarithms_,
  54  *  American Mathematical Society, 1991.
  55  *
  56  *
  57  * ACCURACY:
  58  *
  59  *                      Relative error:
  60  * arithmetic   domain   n   # trials      peak         rms
  61  *    IEEE      0, 1     2     50000      6.2e-16     8.0e-17
  62  *    IEEE      0, 1     3    100000      2.5e-16     6.6e-17
  63  *    IEEE      0, 1     4     30000      1.7e-16     4.9e-17
  64  *    IEEE      0, 1     5     30000      5.1e-16     7.8e-17
  65  *
  66  */
  67
  68 /*
  69 Cephes Math Library Release 2.8:  July, 1999
  70 Copyright 1999 by Stephen L. Moshier
  71 */
  72
  73 #include "mconf.h"
  74 extern double PI;
  75
  76 /* polylog(4, 1-x) = zeta(4) - x zeta(3) + x^2 A4(x)/B4(x)
  77    0 <= x <= 0.125
  78    Theoretical peak absolute error 4.5e-18  */
  79 #if UNK
  80 static double A4[13] = {
  81  3.056144922089490701751E-2,
  82  3.243086484162581557457E-1,
  83  2.877847281461875922565E-1,
  84  7.091267785886180663385E-2,
  85  6.466460072456621248630E-3,
  86  2.450233019296542883275E-4,
  87  4.031655364627704957049E-6,
  88  2.884169163909467997099E-8,
  89  8.680067002466594858347E-11,
  90  1.025983405866370985438E-13,
  91  4.233468313538272640380E-17,
  92  4.959422035066206902317E-21,
  93  1.059365867585275714599E-25,
  94 };
  95 static double B4[12] = {
  96   /* 1.000000000000000000000E0, */
  97  2.821262403600310974875E0,
  98  1.780221124881327022033E0,
  99  3.778888211867875721773E-1,
 100  3.193887040074337940323E-2,
 101  1.161252418498096498304E-3,
 102  1.867362374829870620091E-5,
 103  1.319022779715294371091E-7,
 104  3.942755256555603046095E-10,
 105  4.644326968986396928092E-13,
 106  1.913336021014307074861E-16,
 107  2.240041814626069927477E-20,
 108  4.784036597230791011855E-25,
 109 };
 110 #endif
 111 #if DEC
 112 static short A4[52] = {
 113 0036772,0056001,0016601,0164507,
 114 0037646,0005710,0076603,0176456,
 115 0037623,0054205,0013532,0026476,
 116 0037221,0035252,0101064,0065407,
 117 0036323,0162231,0042033,0107244,
 118 0035200,0073170,0106141,0136543,
 119 0033607,0043647,0163672,0055340,
 120 0031767,0137614,0173376,0072313,
 121 0027676,0160156,0161276,0034203,
 122 0025347,0003752,0123106,0064266,
 123 0022503,0035770,0160173,0177501,
 124 0017273,0056226,0033704,0132530,
 125 0013403,0022244,0175205,0052161,
 126 };
 127 static short B4[48] = {
 128   /*0040200,0000000,0000000,0000000, */
 129 0040464,0107620,0027471,0071672,
 130 0040343,0157111,0025601,0137255,
 131 0037701,0075244,0140412,0160220,
 132 0037002,0151125,0036572,0057163,
 133 0035630,0032452,0050727,0161653,
 134 0034234,0122515,0034323,0172615,
 135 0032415,0120405,0123660,0003160,
 136 0030330,0140530,0161045,0150177,
 137 0026002,0134747,0014542,0002510,
 138 0023134,0113666,0035730,0035732,
 139 0017723,0110343,0041217,0007764,
 140 0014024,0007412,0175575,0160230,
 141 };
 142 #endif
 143 #if IBMPC
 144 static short A4[52] = {
 145 0x3d29,0x23b0,0x4b80,0x3f9f,
 146 0x7fa6,0x0fb0,0xc179,0x3fd4,
 147 0x45a8,0xa2eb,0x6b10,0x3fd2,
 148 0x8d61,0x5046,0x2755,0x3fb2,
 149 0x71d4,0x2883,0x7c93,0x3f7a,
 150 0x37ac,0x118c,0x0ecf,0x3f30,
 151 0x4b5c,0xfcf7,0xe8f4,0x3ed0,
 152 0xce99,0x9edf,0xf7f1,0x3e5e,
 153 0xc710,0xdc57,0xdc0d,0x3dd7,
 154 0xcd17,0x54c8,0xe0fd,0x3d3c,
 155 0x7fe8,0x1c0f,0x677f,0x3c88,
 156 0x96ab,0xc6f8,0x6b92,0x3bb7,
 157 0xaa8e,0x9f50,0x6494,0x3ac0,
 158 };
 159 static short B4[48] = {
 160   /*0x0000,0x0000,0x0000,0x3ff0,*/
 161 0x2e77,0x05e7,0x91f2,0x4006,
 162 0x37d6,0x2570,0x7bc9,0x3ffc,
 163 0x5c12,0x9821,0x2f54,0x3fd8,
 164 0x4bce,0xa7af,0x5a4a,0x3fa0,
 165 0xfc75,0x4a3a,0x06a5,0x3f53,
 166 0x7eb2,0xa71a,0x94a9,0x3ef3,
 167 0x00ce,0xb4f6,0xb420,0x3e81,
 168 0xba10,0x1c44,0x182b,0x3dfb,
 169 0x40a9,0xe32c,0x573c,0x3d60,
 170 0x077b,0xc77b,0x92f6,0x3cab,
 171 0xe1fe,0x6851,0x721c,0x3bda,
 172 0xbc13,0x5f6f,0x81e1,0x3ae2,
 173 };
 174 #endif
 175 #if MIEEE
 176 static short A4[52] = {
 177 0x3f9f,0x4b80,0x23b0,0x3d29,
 178 0x3fd4,0xc179,0x0fb0,0x7fa6,
 179 0x3fd2,0x6b10,0xa2eb,0x45a8,
 180 0x3fb2,0x2755,0x5046,0x8d61,
 181 0x3f7a,0x7c93,0x2883,0x71d4,
 182 0x3f30,0x0ecf,0x118c,0x37ac,
 183 0x3ed0,0xe8f4,0xfcf7,0x4b5c,
 184 0x3e5e,0xf7f1,0x9edf,0xce99,
 185 0x3dd7,0xdc0d,0xdc57,0xc710,
 186 0x3d3c,0xe0fd,0x54c8,0xcd17,
 187 0x3c88,0x677f,0x1c0f,0x7fe8,
 188 0x3bb7,0x6b92,0xc6f8,0x96ab,
 189 0x3ac0,0x6494,0x9f50,0xaa8e,
 190 };
 191 static short B4[48] = {
 192   /*0x3ff0,0x0000,0x0000,0x0000,*/
 193 0x4006,0x91f2,0x05e7,0x2e77,
 194 0x3ffc,0x7bc9,0x2570,0x37d6,
 195 0x3fd8,0x2f54,0x9821,0x5c12,
 196 0x3fa0,0x5a4a,0xa7af,0x4bce,
 197 0x3f53,0x06a5,0x4a3a,0xfc75,
 198 0x3ef3,0x94a9,0xa71a,0x7eb2,
 199 0x3e81,0xb420,0xb4f6,0x00ce,
 200 0x3dfb,0x182b,0x1c44,0xba10,
 201 0x3d60,0x573c,0xe32c,0x40a9,
 202 0x3cab,0x92f6,0xc77b,0x077b,
 203 0x3bda,0x721c,0x6851,0xe1fe,
 204 0x3ae2,0x81e1,0x5f6f,0xbc13,
 205 };
 206 #endif
 207
 208 #ifdef ANSIPROT
 209 extern double spence ( double );
 210 extern double polevl ( double, void *, int );
 211 extern double p1evl ( double, void *, int );
 212 extern double zetac ( double );
 213 extern double pow ( double, double );
 214 extern double powi ( double, int );
 215 extern double log ( double );
 216 extern double fac ( int i );
 217 extern double fabs (double);
 218 double polylog (int, double);
 219 #else
 220 extern double spence(), polevl(), p1evl(), zetac();
 221 extern double pow(), powi(), log();
 222 extern double fac(); /* factorial */
 223 extern double fabs();
 224 double polylog();
 225 #endif
 226 extern double MACHEP;
 227
 228 double
 229 polylog (n, x)
 230      int n;
 231      double x;
 232 {
 233   double h, k, p, s, t, u, xc, z;
 234   int i, j;
 235
 236 /*  This recurrence provides formulas for n < 2.
 237
 238     d                 1
 239     --   Li (x)  =   ---  Li   (x)  .
 240     dx     n          x     n-1
 241
 242 */
 243
 244   if (n == -1)
 245     {
 246       p  = 1.0 - x;
 247       u = x / p;
 248       s = u * u + u;
 249       return s;
 250     }
 251
 252   if (n == 0)
 253     {
 254       s = x / (1.0 - x);
 255       return s;
 256     }
 257
 258   /* Not implemented for n < -1.
 259      Not defined for x > 1.  Use cpolylog if you need that.  */
 260   if (x > 1.0 || n < -1)
 261     {
 262       mtherr("polylog", DOMAIN);
 263       return 0.0;
 264     }
 265
 266   if (n == 1)
 267     {
 268       s = -log (1.0 - x);
 269       return s;
 270     }
 271
 272   /* Argument +1 */
 273   if (x == 1.0 && n > 1)
 274     {
 275       s = zetac ((double) n) + 1.0;
 276       return s;
 277     }
 278
 279   /* Argument -1.
 280                         1-n
 281      Li (-z)  = - (1 - 2   ) Li (z)
 282        n                       n
 283    */
 284   if (x == -1.0 && n > 1)
 285     {
 286       /* Li_n(1) = zeta(n) */
 287       s = zetac ((double) n) + 1.0;
 288       s = s * (powi (2.0, 1 - n) - 1.0);
 289       return s;
 290     }
 291
 292 /*  Inversion formula:
 293  *                                                   [n/2]   n-2r
 294  *                n                  1     n           -  log    (z)
 295  *  Li (-z) + (-1)  Li (-1/z)  =  - --- log (z)  +  2  >  ----------- Li  (-1)
 296  *    n               n              n!                -   (n - 2r)!    2r
 297  *                                                    r=1
 298  */
 299   if (x < -1.0 && n > 1)
 300     {
 301       double q, w;
 302       int r;
 303
 304       w = log (-x);
 305       s = 0.0;
 306       for (r = 1; r <= n / 2; r++)
 307         {
 308           j = 2 * r;
 309           p = polylog (j, -1.0);
 310           j = n - j;
 311           if (j == 0)
 312             {
 313               s = s + p;
 314               break;
 315             }
 316           q = (double) j;
 317           q = pow (w, q) * p / fac (j);
 318           s = s + q;
 319         }
 320       s = 2.0 * s;
 321       q = polylog (n, 1.0 / x);
 322       if (n & 1)
 323         q = -q;
 324       s = s - q;
 325       s = s - pow (w, (double) n) / fac (n);
 326       return s;
 327     }
 328
 329   if (n == 2)
 330     {
 331       if (x < 0.0 || x > 1.0)
 332         return (spence (1.0 - x));
 333     }
 334
 335
 336
 337   /*  The power series converges slowly when x is near 1.  For n = 3, this
 338       identity helps:
 339
 340       Li (-x/(1-x)) + Li (1-x) + Li (x)
 341         3               3          3
 342                      2                               2                 3
 343        = Li (1) + (pi /6) log(1-x) - (1/2) log(x) log (1-x) + (1/6) log (1-x)
 344            3
 345   */
 346
 347   if (n == 3)
 348     {
 349       p = x * x * x;
 350       if (x > 0.8)
 351         {
 352           u = log(x);
 353           s = p / 6.0;
 354           xc = 1.0 - x;
 355           s = s - 0.5 * u * u * log(xc);
 356           s = s + PI * PI * u / 6.0;
 357           s = s - polylog (3, -xc/x);
 358           s = s - polylog (3, xc);
 359           s = s + zetac(3.0);
 360           s = s + 1.0;
 361           return s;
 362         }
 363       /* Power series  */
 364       t = p / 27.0;
 365       t = t + .125 * x * x;
 366       t = t + x;
 367
 368       s = 0.0;
 369       k = 4.0;
 370       do
 371         {
 372           p = p * x;
 373           h = p / (k * k * k);
 374           s = s + h;
 375           k += 1.0;
 376         }
 377       while (fabs(h/s) > 1.1e-16);
 378       return (s + t);
 379     }
 380
 381 if (n == 4)
 382   {
 383     if (x >= 0.875)
 384       {
 385         u = 1.0 - x;
 386         s = polevl(u, A4, 12) / p1evl(u, B4, 12);
 387         s =  s * u * u - 1.202056903159594285400 * u;
 388         s +=  1.0823232337111381915160;
 389         return s;
 390       }
 391     goto pseries;
 392   }
 393
 394
 395   if (x < 0.75)
 396     goto pseries;
 397
 398
 399 /*  This expansion in powers of log(x) is especially useful when
 400     x is near 1.
 401
 402     See also the pari gp calculator.
 403
 404                       inf                  j
 405                        -    z(n-j) (log(x))
 406     polylog(n,x)  =    >   -----------------
 407                        -           j!
 408                       j=0
 409
 410       where
 411
 412       z(j) = Riemann zeta function (j), j != 1
 413
 414                               n-1
 415                                -
 416       z(1) =  -log(-log(x)) +  >  1/k
 417                                -
 418                               k=1
 419   */
 420
 421   z = log(x);
 422   h = -log(-z);
 423   for (i = 1; i < n; i++)
 424     h = h + 1.0/i;
 425   p = 1.0;
 426   s = zetac((double)n) + 1.0;
 427   for (j=1; j<=n+1; j++)
 428   {
 429     p = p * z / j;
 430     if (j == n-1)
 431       s = s + h * p;
 432     else
 433       s = s + (zetac((double)(n-j)) + 1.0) * p;
 434   }
 435   j = n + 3;
 436   z = z * z;
 437   for(;;)
 438     {
 439       p = p * z / ((j-1)*j);
 440       h = (zetac((double)(n-j)) + 1.0);
 441       h = h * p;
 442       s = s + h;
 443       if (fabs(h/s) < MACHEP)
 444         break;
 445       j += 2;
 446     }
 447   return s;
 448
 449
 450 pseries:
 451
 452   p = x * x * x;
 453   k = 3.0;
 454   s = 0.0;
 455   do
 456     {
 457       p = p * x;
 458       k += 1.0;
 459       h = p / powi(k, n);
 460       s = s + h;
 461     }
 462   while (fabs(h/s) > MACHEP);
 463   s += x * x * x / powi(3.0,n);
 464   s += x * x / powi(2.0,n);
 465   s += x;
 466   return s;
 467 }