From 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 22 Nov 2001 14:04:29 +0000 Subject: Totally rework the math library, this time based on the MacOs X math library (which is itself based on the math lib from FreeBSD). -Erik --- libm/double/psi.c | 201 ------------------------------------------------------ 1 file changed, 201 deletions(-) delete mode 100644 libm/double/psi.c (limited to 'libm/double/psi.c') diff --git a/libm/double/psi.c b/libm/double/psi.c deleted file mode 100644 index 6da2aa0c2..000000000 --- a/libm/double/psi.c +++ /dev/null @@ -1,201 +0,0 @@ -/* psi.c - * - * Psi (digamma) function - * - * - * SYNOPSIS: - * - * double x, y, psi(); - * - * y = psi( 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: - * Relative error (except absolute when |psi| < 1): - * arithmetic domain # trials peak rms - * DEC 0,30 2500 1.7e-16 2.0e-17 - * IEEE 0,30 30000 1.3e-15 1.4e-16 - * IEEE -30,0 40000 1.5e-15 2.2e-16 - * - * ERROR MESSAGES: - * message condition value returned - * psi singularity x integer <=0 MAXNUM - */ - -/* -Cephes Math Library Release 2.8: June, 2000 -Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier -*/ - -#include - -#ifdef UNK -static double A[] = { - 8.33333333333333333333E-2, --2.10927960927960927961E-2, - 7.57575757575757575758E-3, --4.16666666666666666667E-3, - 3.96825396825396825397E-3, --8.33333333333333333333E-3, - 8.33333333333333333333E-2 -}; -#endif - -#ifdef DEC -static unsigned short A[] = { -0037252,0125252,0125252,0125253, -0136654,0145314,0126312,0146255, -0036370,0037017,0101740,0174076, -0136210,0104210,0104210,0104211, -0036202,0004040,0101010,0020202, -0136410,0104210,0104210,0104211, -0037252,0125252,0125252,0125253 -}; -#endif - -#ifdef IBMPC -static unsigned short A[] = { -0x5555,0x5555,0x5555,0x3fb5, -0x5996,0x9599,0x9959,0xbf95, -0x1f08,0xf07c,0x07c1,0x3f7f, -0x1111,0x1111,0x1111,0xbf71, -0x0410,0x1041,0x4104,0x3f70, -0x1111,0x1111,0x1111,0xbf81, -0x5555,0x5555,0x5555,0x3fb5 -}; -#endif - -#ifdef MIEEE -static unsigned short A[] = { -0x3fb5,0x5555,0x5555,0x5555, -0xbf95,0x9959,0x9599,0x5996, -0x3f7f,0x07c1,0xf07c,0x1f08, -0xbf71,0x1111,0x1111,0x1111, -0x3f70,0x4104,0x1041,0x0410, -0xbf81,0x1111,0x1111,0x1111, -0x3fb5,0x5555,0x5555,0x5555 -}; -#endif - -#define EUL 0.57721566490153286061 - -#ifdef ANSIPROT -extern double floor ( double ); -extern double log ( double ); -extern double tan ( double ); -extern double polevl ( double, void *, int ); -#else -double floor(), log(), tan(), polevl(); -#endif -extern double PI, MAXNUM; - - -double psi(x) -double x; -{ -double p, q, nz, s, w, y, z; -int i, n, negative; - -negative = 0; -nz = 0.0; - -if( x <= 0.0 ) - { - negative = 1; - q = x; - p = floor(q); - if( p == q ) - { - mtherr( "psi", SING ); - return( MAXNUM ); - } -/* Remove the zeros of tan(PI x) - * by subtracting the nearest integer from x - */ - nz = q - p; - if( nz != 0.5 ) - { - if( nz > 0.5 ) - { - p += 1.0; - nz = q - p; - } - nz = PI/tan(PI*nz); - } - else - { - nz = 0.0; - } - x = 1.0 - x; - } - -/* check for positive integer up to 10 */ -if( (x <= 10.0) && (x == floor(x)) ) - { - y = 0.0; - n = x; - for( i=1; i