summaryrefslogtreecommitdiff
path: root/libm/double/psi.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/double/psi.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
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
Diffstat (limited to 'libm/double/psi.c')
-rw-r--r--libm/double/psi.c201
1 files changed, 0 insertions, 201 deletions
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 <math.h>
-
-#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<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.0e17 )
- {
- z = 1.0/(s * s);
- y = z * polevl( z, A, 6 );
- }
-else
- y = 0.0;
-
-y = log(s) - (0.5/s) - y - w;
-
-done:
-
-if( negative )
- {
- y -= nz;
- }
-
-return(y);
-}