view release on metacpan or  search on metacpan

libmd/polylog.c  view on Meta::CPAN

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     n-1

*/

  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.
                        1-n
     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]   n-2r
 *                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/(1-x)) + Li (1-x) + Li (x)
        3               3          3
                     2                               2                 3
       = Li (1) + (pi /6) md_log(1-x) - (1/2) md_log(x) md_log (1-x) + (1/6) md_log (1-x)
           3
  */

  if (n == 3)
    {
      p = x * x * x;
      if (x > 0.8)
	{
	  u = md_log(x);