diff options
Diffstat (limited to 'libm/float/psif.c')
-rw-r--r-- | libm/float/psif.c | 153 |
1 files changed, 0 insertions, 153 deletions
diff --git a/libm/float/psif.c b/libm/float/psif.c deleted file mode 100644 index 2d9187c67..000000000 --- a/libm/float/psif.c +++ /dev/null @@ -1,153 +0,0 @@ -/* 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); -} |