// psi.cpp -- psi function for real arguments.
//      Algorithms and coefficient values from "Computation of Special
//      Functions", Zhang and Jin, John Wiley and Sons, 1996.
//
//  (C) 2003, C. Bond. All rights reserved.
//
//  Returns psi function for real argument 'x'.
//  NOTE: Returns 1e308 for argument 0 or a negative integer.
//
//#include <iostream.h>
#include <math.h>
#define el 0.5772156649015329

double psi ( double x ) {
	double s, ps, xa, x2;
	int n, k;
	static double a[] = {
		-0.8333333333333e-01,
		0.83333333333333333e-02,
		-0.39682539682539683e-02,
		0.41666666666666667e-02,
		-0.75757575757575758e-02,
		0.21092796092796093e-01,
		-0.83333333333333333e-01,
		0.4432598039215686
	};

	xa = fabs ( x );
	s = 0.0;
	if ( ( x == ( int ) x ) && ( x <= 0.0 ) ) {
		ps = 1e308;
		return ps;
	}
	if ( xa == ( int ) xa ) {
		n = xa;
		for ( k = 1; k < n; k++ ) {
			s += 1.0 / k;
		}
		ps =  s - el;
	} else if ( ( xa + 0.5 ) == ( ( int ) ( xa + 0.5 ) ) ) {
		n = xa - 0.5;
		for ( k = 1; k <= n; k++ ) {
			s += 1.0 / ( 2.0 * k - 1.0 );
		}
		ps = 2.0 * s - el - 1.386294361119891;
	} else {
		if ( xa < 10.0 ) {
			n = 10 - ( int ) xa;
			for ( k = 0; k < n; k++ ) {
				s += 1.0 / ( xa + k );
			}
			xa += n;
		}
		x2 = 1.0 / ( xa * xa );
		ps = log ( xa ) - 0.5 / xa + x2 * ( ( ( ( ( ( ( a[7] * x2 + a[6] ) * x2 + a[5] ) * x2 +
		                                        a[4] ) * x2 + a[3] ) * x2 + a[2] ) * x2 + a[1] ) * x2 + a[0] );
		ps -= s;
	}
	if ( x < 0.0 )
		ps = ps - M_PI * cos ( M_PI * x ) / sin ( M_PI * x ) - 1.0 / x;
	return ps;
}