// 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;
}