From 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 Mon Sep 17 00:00:00 2001
From: Eric Andersen <andersen@codepoet.org>
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/float/psif.c | 153 ------------------------------------------------------
 1 file changed, 153 deletions(-)
 delete mode 100644 libm/float/psif.c

(limited to 'libm/float/psif.c')

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);
-}
-- 
cgit v1.2.3