/*							psif.c
 *
 *	Psi (digamma) function
 *
 *
 * SYNOPSIS:
 *
 * float x, y, psif();
 *
 * y = psif( x );
 *
 *
 * DESCRIPTION:
 *
 *              d      -
 *   psi(x)  =  -- ln | (x)
 *              dx
 *
 * is the logarithmic derivative of the gamma function.
 * For integer x,
 *                   n-1
 *                    -
 * psi(n) = -EUL  +   >  1/k.
 *                    -
 *                   k=1
 *
 * This formula is used for 0 < n <= 10.  If x is negative, it
 * is transformed to a positive argument by the reflection
 * formula  psi(1-x) = psi(x) + pi cot(pi x).
 * For general positive x, the argument is made greater than 10
 * using the recurrence  psi(x+1) = psi(x) + 1/x.
 * Then the following asymptotic expansion is applied:
 *
 *                           inf.   B
 *                            -      2k
 * psi(x) = log(x) - 1/2x -   >   -------
 *                            -        2k
 *                           k=1   2k x
 *
 * where the B2k are Bernoulli numbers.
 *
 * ACCURACY:
 *    Absolute error,  relative when |psi| > 1 :
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -33,0        30000      8.2e-7      1.2e-7
 *    IEEE      0,33        100000      7.3e-7      7.7e-8
 *
 * ERROR MESSAGES:
 *     message         condition      value returned
 * psi singularity    x integer <=0      MAXNUMF
 */

/*
Cephes Math Library Release 2.2:  June, 1992
Copyright 1984, 1987, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

#include <math.h>


static float A[] = {
-4.16666666666666666667E-3,
 3.96825396825396825397E-3,
-8.33333333333333333333E-3,
 8.33333333333333333333E-2
};


#define EUL 0.57721566490153286061

extern float PIF, MAXNUMF;



float floorf(float), logf(float), tanf(float);
float polevlf(float, float *, int);

float psif(float xx)
{
float p, q, nz, x, s, w, y, z;
int i, n, negative;


x = xx;
nz = 0.0;
negative = 0;
if( x <= 0.0 )
	{
	negative = 1;
	q = x;
	p = floorf(q);
	if( p == q )
		{
		mtherr( "psif", SING );
		return( MAXNUMF );
		}
	nz = q - p;
	if( nz != 0.5 )
		{
		if( nz > 0.5 )
			{
			p += 1.0;
			nz = q - p;
			}
		nz = PIF/tanf(PIF*nz);
		}
	else
		{
		nz = 0.0;
		}
	x = 1.0 - x;
	}

/* check for positive integer up to 10 */
if( (x <= 10.0) && (x == floorf(x)) )
	{
	y = 0.0;
	n = x;
	for( i=1; i<n; i++ )
		{
		w = i;
		y += 1.0/w;
		}
	y -= EUL;
	goto done;
	}

s = x;
w = 0.0;
while( s < 10.0 )
	{
	w += 1.0/s;
	s += 1.0;
	}

if( s < 1.0e8 )
	{
	z = 1.0/(s * s);
	y = z * polevlf( z, A, 3 );
	}
else
	y = 0.0;

y = logf(s)  -  (0.5/s)  -  y  -  w;

done:
if( negative )
	{
	y -= nz;
	}
return(y);
}