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