lib/libm/psi.c

   1 /*                                                      psi.c
   2  *
   3  *      Psi (digamma) function
   4  *
   5  *
   6  * SYNOPSIS:
   7  *
   8  * double x, y, psi();
   9  *
  10  * y = psi( x );
  11  *
  12  *
  13  * DESCRIPTION:
  14  *
  15  *              d      -
  16  *   psi(x)  =  -- ln | (x)
  17  *              dx
  18  *
  19  * is the logarithmic derivative of the gamma function.
  20  * For integer x,
  21  *                   n-1
  22  *                    -
  23  * psi(n) = -EUL  +   >  1/k.
  24  *                    -
  25  *                   k=1
  26  *
  27  * This formula is used for 0 < n <= 10.  If x is negative, it
  28  * is transformed to a positive argument by the reflection
  29  * formula  psi(1-x) = psi(x) + pi cot(pi x).
  30  * For general positive x, the argument is made greater than 10
  31  * using the recurrence  psi(x+1) = psi(x) + 1/x.
  32  * Then the following asymptotic expansion is applied:
  33  *
  34  *                           inf.   B
  35  *                            -      2k
  36  * psi(x) = log(x) - 1/2x -   >   -------
  37  *                            -        2k
  38  *                           k=1   2k x
  39  *
  40  * where the B2k are Bernoulli numbers.
  41  *
  42  * ACCURACY:
  43  *    Relative error (except absolute when |psi| < 1):
  44  * arithmetic   domain     # trials      peak         rms
  45  *    DEC       0,30         2500       1.7e-16     2.0e-17
  46  *    IEEE      0,30        30000       1.3e-15     1.4e-16
  47  *    IEEE      -30,0       40000       1.5e-15     2.2e-16
  48  *
  49  * ERROR MESSAGES:
  50  *     message         condition      value returned
  51  * psi singularity    x integer <=0      MAXNUM
  52  */
  53 \f
  54 /*
  55 Cephes Math Library Release 2.8:  June, 2000
  56 Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier
  57 */
  58
  59 #include "mconf.h"
  60
  61 #ifdef UNK
  62 static double A[] = {
  63  8.33333333333333333333E-2,
  64 -2.10927960927960927961E-2,
  65  7.57575757575757575758E-3,
  66 -4.16666666666666666667E-3,
  67  3.96825396825396825397E-3,
  68 -8.33333333333333333333E-3,
  69  8.33333333333333333333E-2
  70 };
  71 #endif
  72
  73 #ifdef DEC
  74 static unsigned short A[] = {
  75 0037252,0125252,0125252,0125253,
  76 0136654,0145314,0126312,0146255,
  77 0036370,0037017,0101740,0174076,
  78 0136210,0104210,0104210,0104211,
  79 0036202,0004040,0101010,0020202,
  80 0136410,0104210,0104210,0104211,
  81 0037252,0125252,0125252,0125253
  82 };
  83 #endif
  84
  85 #ifdef IBMPC
  86 static unsigned short A[] = {
  87 0x5555,0x5555,0x5555,0x3fb5,
  88 0x5996,0x9599,0x9959,0xbf95,
  89 0x1f08,0xf07c,0x07c1,0x3f7f,
  90 0x1111,0x1111,0x1111,0xbf71,
  91 0x0410,0x1041,0x4104,0x3f70,
  92 0x1111,0x1111,0x1111,0xbf81,
  93 0x5555,0x5555,0x5555,0x3fb5
  94 };
  95 #endif
  96
  97 #ifdef MIEEE
  98 static unsigned short A[] = {
  99 0x3fb5,0x5555,0x5555,0x5555,
 100 0xbf95,0x9959,0x9599,0x5996,
 101 0x3f7f,0x07c1,0xf07c,0x1f08,
 102 0xbf71,0x1111,0x1111,0x1111,
 103 0x3f70,0x4104,0x1041,0x0410,
 104 0xbf81,0x1111,0x1111,0x1111,
 105 0x3fb5,0x5555,0x5555,0x5555
 106 };
 107 #endif
 108
 109 #define EUL 0.57721566490153286061
 110
 111 #ifdef ANSIPROT
 112 extern double floor ( double );
 113 extern double log ( double );
 114 extern double tan ( double );
 115 extern double polevl ( double, void *, int );
 116 #else
 117 double floor(), log(), tan(), polevl();
 118 #endif
 119 extern double PI, MAXNUM;
 120
 121
 122 double psi(x)
 123 double x;
 124 {
 125 double p, q, nz, s, w, y, z;
 126 int i, n, negative;
 127
 128 negative = 0;
 129 nz = 0.0;
 130
 131 if( x <= 0.0 )
 132         {
 133         negative = 1;
 134         q = x;
 135         p = floor(q);
 136         if( p == q )
 137                 {
 138                 mtherr( "psi", SING );
 139                 return( MAXNUM );
 140                 }
 141 /* Remove the zeros of tan(PI x)
 142  * by subtracting the nearest integer from x
 143  */
 144         nz = q - p;
 145         if( nz != 0.5 )
 146                 {
 147                 if( nz > 0.5 )
 148                         {
 149                         p += 1.0;
 150                         nz = q - p;
 151                         }
 152                 nz = PI/tan(PI*nz);
 153                 }
 154         else
 155                 {
 156                 nz = 0.0;
 157                 }
 158         x = 1.0 - x;
 159         }
 160
 161 /* check for positive integer up to 10 */
 162 if( (x <= 10.0) && (x == floor(x)) )
 163         {
 164         y = 0.0;
 165         n = x;
 166         for( i=1; i<n; i++ )
 167                 {
 168                 w = i;
 169                 y += 1.0/w;
 170                 }
 171         y -= EUL;
 172         goto done;
 173         }
 174
 175 s = x;
 176 w = 0.0;
 177 while( s < 10.0 )
 178         {
 179         w += 1.0/s;
 180         s += 1.0;
 181         }
 182
 183 if( s < 1.0e17 )
 184         {
 185         z = 1.0/(s * s);
 186         y = z * polevl( z, A, 6 );
 187         }
 188 else
 189         y = 0.0;
 190
 191 y = log(s)  -  (0.5/s)  -  y  -  w;
 192
 193 done:
 194
 195 if( negative )
 196         {
 197         y -= nz;
 198         }
 199
 200 return(y);
 201 }