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

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