uclibc-ng/libm/float/spencef.c

/*							spencef.c
 *
 *	Dilogarithm
 *
 *
 *
 * SYNOPSIS:
 *
 * float x, y, spencef();
 *
 * y = spencef( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Computes the integral
 *
 *                    x
 *                    -
 *                   | | log t
 * spence(x)  =  -   |   ----- dt
 *                 | |   t - 1
 *                  -
 *                  1
 *
 * for x >= 0.  A rational approximation gives the integral in
 * the interval (0.5, 1.5).  Transformation formulas for 1/x
 * and 1-x are employed outside the basic expansion range.
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      0,4         30000       4.4e-7      6.3e-8
 *
 *
 */

/*							spence.c */


/*
Cephes Math Library Release 2.1:  January, 1989
Copyright 1985, 1987, 1989 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

#include <math.h>

static float A[8] = {
  4.65128586073990045278E-5,
  7.31589045238094711071E-3,
  1.33847639578309018650E-1,
  8.79691311754530315341E-1,
  2.71149851196553469920E0,
  4.25697156008121755724E0,
  3.29771340985225106936E0,
  1.00000000000000000126E0,
};
static float B[8] = {
  6.90990488912553276999E-4,
  2.54043763932544379113E-2,
  2.82974860602568089943E-1,
  1.41172597751831069617E0,
  3.63800533345137075418E0,
  5.03278880143316990390E0,
  3.54771340985225096217E0,
  9.99999999999999998740E-1,
};

extern float PIF, MACHEPF;

/* pi * pi / 6 */
#define PIFS 1.64493406684822643647


float logf(float), polevlf(float, float *, int);
float spencef(float xx)
{
float x, w, y, z;
int flag;

x = xx;
if( x < 0.0 )
	{
	mtherr( "spencef", DOMAIN );
	return(0.0);
	}

if( x == 1.0 )
	return( 0.0 );

if( x == 0.0 )
	return( PIFS );

flag = 0;

if( x > 2.0 )
	{
	x = 1.0/x;
	flag |= 2;
	}

if( x > 1.5 )
	{
	w = (1.0/x) - 1.0;
	flag |= 2;
	}

else if( x < 0.5 )
	{
	w = -x;
	flag |= 1;
	}

else
	w = x - 1.0;


y = -w * polevlf( w, A, 7) / polevlf( w, B, 7 );

if( flag & 1 )
	y = PIFS - logf(x) * logf(1.0-x) - y;

if( flag & 2 )
	{
	z = logf(x);
	y = -0.5 * z * z  -  y;
	}

return( y );
}