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
54 /*
55 Cephes Math Library Release 2.2: July, 1992
56 Copyright 1984, 1987, 1992 by Stephen L. Moshier
57 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
58 */
59
60 #include "mconf.h"
61 #include "cephes.h"
62
63 #ifdef UNK
64 static double A[] = {
65 8.33333333333333333333E-2,
66 -2.10927960927960927961E-2,
67 7.57575757575757575758E-3,
68 -4.16666666666666666667E-3,
69 3.96825396825396825397E-3,
70 -8.33333333333333333333E-3,
71 8.33333333333333333333E-2
72 };
73 #endif
74
75 #ifdef DEC
76 static unsigned short A[] = {
77 0037252,0125252,0125252,0125253,
78 0136654,0145314,0126312,0146255,
79 0036370,0037017,0101740,0174076,
80 0136210,0104210,0104210,0104211,
81 0036202,0004040,0101010,0020202,
82 0136410,0104210,0104210,0104211,
83 0037252,0125252,0125252,0125253
84 };
85 #endif
86
87 #ifdef IBMPC
88 static unsigned short A[] = {
89 0x5555,0x5555,0x5555,0x3fb5,
90 0x5996,0x9599,0x9959,0xbf95,
91 0x1f08,0xf07c,0x07c1,0x3f7f,
92 0x1111,0x1111,0x1111,0xbf71,
93 0x0410,0x1041,0x4104,0x3f70,
94 0x1111,0x1111,0x1111,0xbf81,
95 0x5555,0x5555,0x5555,0x3fb5
96 };
97 #endif
98
99 #ifdef MIEEE
100 static unsigned short A[] = {
101 0x3fb5,0x5555,0x5555,0x5555,
102 0xbf95,0x9959,0x9599,0x5996,
103 0x3f7f,0x07c1,0xf07c,0x1f08,
104 0xbf71,0x1111,0x1111,0x1111,
105 0x3f70,0x4104,0x1041,0x0410,
106 0xbf81,0x1111,0x1111,0x1111,
107 0x3fb5,0x5555,0x5555,0x5555
108 };
109 #endif
110
111 #define EUL 0.57721566490153286061
112
113 extern double PI, MAXNUM;
114
115
psi(x)116 double psi(x)
117 double x;
118 {
119 double p, q, nz, s, w, y, z;
120 int i, n, negative;
121
122 negative = 0;
123 nz = 0.0;
124
125 if( x <= 0.0 )
126 {
127 negative = 1;
128 q = x;
129 p = floor(q);
130 if( p == q )
131 {
132 mtherr( "psi", SING );
133 return( MAXNUM );
134 }
135 /* Remove the zeros of tan(PI x)
136 * by subtracting the nearest integer from x
137 */
138 nz = q - p;
139 if( nz != 0.5 )
140 {
141 if( nz > 0.5 )
142 {
143 p += 1.0;
144 nz = q - p;
145 }
146 nz = PI/tan(PI*nz);
147 }
148 else
149 {
150 nz = 0.0;
151 }
152 x = 1.0 - x;
153 }
154
155 /* check for positive integer up to 10 */
156 if( (x <= 10.0) && (x == floor(x)) )
157 {
158 y = 0.0;
159 n = x;
160 for( i=1; i<n; i++ )
161 {
162 w = i;
163 y += 1.0/w;
164 }
165 y -= EUL;
166 goto done;
167 }
168
169 s = x;
170 w = 0.0;
171 while( s < 10.0 )
172 {
173 w += 1.0/s;
174 s += 1.0;
175 }
176
177 if( s < 1.0e17 )
178 {
179 z = 1.0/(s * s);
180 y = z * polevl( z, A, 6 );
181 }
182 else
183 y = 0.0;
184
185 y = log(s) - (0.5/s) - y - w;
186
187 done:
188
189 if( negative )
190 {
191 y -= nz;
192 }
193
194 return(y);
195 }
196