root/library/bdm/math/psi.cpp @ 1187

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