1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 1998 Ross Ihaka
4 * Copyright (C) 2000-2007 the R Core Team
5 * Copyright (C) 2004-2009 The R Foundation
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, a copy is available at
19 * http://www.r-project.org/Licenses/
20 *
21 * SYNOPSIS
22 *
23 * #include <Rmath.h>
24 * void dpsifn(double x, int n, int kode, int m,
25 * double *ans, int *nz, int *ierr)
26 * double digamma(double x);
27 * double trigamma(double x)
28 * double tetragamma(double x)
29 * double pentagamma(double x)
30 *
31 * DESCRIPTION
32 *
33 * Compute the derivatives of the psi function
34 * and polygamma functions.
35 *
36 * The following definitions are used in dpsifn:
37 *
38 * Definition 1
39 *
40 * psi(x) = d/dx (ln(gamma(x)), the first derivative of
41 * the log gamma function.
42 *
43 * Definition 2
44 * k k
45 * psi(k,x) = d /dx (psi(x)), the k-th derivative
46 * of psi(x).
47 *
48 *
49 * "dpsifn" computes a sequence of scaled derivatives of
50 * the psi function; i.e. for fixed x and m it computes
51 * the m-member sequence
52 *
53 * (-1)^(k+1) / gamma(k+1) * psi(k,x)
54 * for k = n,...,n+m-1
55 *
56 * where psi(k,x) is as defined above. For kode=1, dpsifn
57 * returns the scaled derivatives as described. kode=2 is
58 * operative only when k=0 and in that case dpsifn returns
59 * -psi(x) + ln(x). That is, the logarithmic behavior for
60 * large x is removed when kode=2 and k=0. When sums or
61 * differences of psi functions are computed the logarithmic
62 * terms can be combined analytically and computed separately
63 * to help retain significant digits.
64 *
65 * Note that dpsifn(x, 0, 1, 1, ans) results in ans = -psi(x).
66 *
67 * INPUT
68 *
69 * x - argument, x > 0.
70 *
71 * n - first member of the sequence, 0 <= n <= 100
72 * n == 0 gives ans(1) = -psi(x) for kode=1
73 * -psi(x)+ln(x) for kode=2
74 *
75 * kode - selection parameter
76 * kode == 1 returns scaled derivatives of the
77 * psi function.
78 * kode == 2 returns scaled derivatives of the
79 * psi function except when n=0. In this case,
80 * ans(1) = -psi(x) + ln(x) is returned.
81 *
82 * m - number of members of the sequence, m >= 1
83 *
84 * OUTPUT
85 *
86 * ans - a vector of length at least m whose first m
87 * components contain the sequence of derivatives
88 * scaled according to kode.
89 *
90 * nz - underflow flag
91 * nz == 0, a normal return
92 * nz != 0, underflow, last nz components of ans are
93 * set to zero, ans(m-k+1)=0.0, k=1,...,nz
94 *
95 * ierr - error flag
96 * ierr=0, a normal return, computation completed
97 * ierr=1, input error, no computation
98 * ierr=2, overflow, x too small or n+m-1 too
99 * large or both
100 * ierr=3, error, n too large. dimensioned
101 * array trmr(nmax) is not large enough for n
102 *
103 * The nominal computational accuracy is the maximum of unit
104 * roundoff (d1mach(4)) and 1e-18 since critical constants
105 * are given to only 18 digits.
106 *
107 * The basic method of evaluation is the asymptotic expansion
108 * for large x >= xmin followed by backward recursion on a two
109 * term recursion relation
110 *
111 * w(x+1) + x^(-n-1) = w(x).
112 *
113 * this is supplemented by a series
114 *
115 * sum( (x+k)^(-n-1) , k=0,1,2,... )
116 *
117 * which converges rapidly for large n. both xmin and the
118 * number of terms of the series are calculated from the unit
119 * roundoff of the machine environment.
120 *
121 * AUTHOR
122 *
123 * Amos, D. E. (Fortran)
124 * Ross Ihaka (C Translation)
125 * Martin Maechler (x < 0, and psigamma())
126 *
127 * REFERENCES
128 *
129 * Handbook of Mathematical Functions,
130 * National Bureau of Standards Applied Mathematics Series 55,
131 * Edited by M. Abramowitz and I. A. Stegun, equations 6.3.5,
132 * 6.3.18, 6.4.6, 6.4.9 and 6.4.10, pp.258-260, 1964.
133 *
134 * D. E. Amos, (1983). "A Portable Fortran Subroutine for
135 * Derivatives of the Psi Function", Algorithm 610,
136 * TOMS 9(4), pp. 494-502.
137 *
138 * Routines called: Rf_d1mach, Rf_i1mach.
139 */
140
141 #include "nmath.h"
142 #ifdef MATHLIB_STANDALONE
143 #include <errno.h>
144 #endif
145
146 #define n_max (100)
147
148 /* From R, currently only used for kode = 1, m = 1, n in {0,1,2,3} : */
dpsifn(double x,int n,int kode,int m,double * ans,int * nz,int * ierr)149 void dpsifn(double x, int n, int kode, int m, double *ans, int *nz, int *ierr)
150 {
151 const static double bvalues[] = { /* Bernoulli Numbers */
152 1.00000000000000000e+00,
153 -5.00000000000000000e-01,
154 1.66666666666666667e-01,
155 -3.33333333333333333e-02,
156 2.38095238095238095e-02,
157 -3.33333333333333333e-02,
158 7.57575757575757576e-02,
159 -2.53113553113553114e-01,
160 1.16666666666666667e+00,
161 -7.09215686274509804e+00,
162 5.49711779448621554e+01,
163 -5.29124242424242424e+02,
164 6.19212318840579710e+03,
165 -8.65802531135531136e+04,
166 1.42551716666666667e+06,
167 -2.72982310678160920e+07,
168 6.01580873900642368e+08,
169 -1.51163157670921569e+10,
170 4.29614643061166667e+11,
171 -1.37116552050883328e+13,
172 4.88332318973593167e+14,
173 -1.92965793419400681e+16
174 };
175
176 int i, j, k, mm, mx, nn, np, nx, fn;
177 double arg, den, elim, eps, fln, fx, rln, rxsq,
178 r1m4, r1m5, s, slope, t, ta, tk, tol, tols, tss, tst,
179 tt, t1, t2, wdtol, xdmln, xdmy, xinc, xln = 0.0 /* -Wall */,
180 xm, xmin, xq, yint;
181 double trm[23], trmr[n_max + 1];
182
183 *ierr = 0;
184 if (n < 0 || kode < 1 || kode > 2 || m < 1) {
185 *ierr = 1;
186 return;
187 }
188 if (x <= 0.) {
189 /* use Abramowitz & Stegun 6.4.7 "Reflection Formula"
190 * psi(k, x) = (-1)^k psi(k, 1-x) - pi^{n+1} (d/dx)^n cot(x)
191 */
192 if (x == (long)x) {
193 /* non-positive integer : +Inf or NaN depends on n */
194 for(j=0; j < m; j++) /* k = j + n : */
195 ans[j] = ((j+n) % 2) ? ML_POSINF : ML_NAN;
196 return;
197 }
198 /* This could cancel badly */
199 dpsifn(1. - x, n, /*kode = */ 1, m, ans, nz, ierr);
200 /* ans[j] == (-1)^(k+1) / gamma(k+1) * psi(k, 1 - x)
201 * for j = 0:(m-1) , k = n + j
202 */
203
204 /* Cheat for now: only work for m = 1, n in {0,1,2,3} : */
205 if(m > 1 || n > 3) {/* doesn't happen for digamma() .. pentagamma() */
206 /* not yet implemented */
207 *ierr = 4; return;
208 }
209 x *= M_PI; /* pi * x */
210 if (n == 0)
211 tt = cos(x)/sin(x);
212 else if (n == 1)
213 tt = -1/pow(sin(x),2);
214 else if (n == 2)
215 tt = 2*cos(x)/pow(sin(x),3);
216 else if (n == 3)
217 tt = -2*(2*pow(cos(x),2) + 1)/pow(sin(x),4);
218 else /* can not happen! */
219 tt = ML_NAN;
220 /* end cheat */
221
222 s = (n % 2) ? -1. : 1.;/* s = (-1)^n */
223 /* t := pi^(n+1) * d_n(x) / gamma(n+1) , where
224 * d_n(x) := (d/dx)^n cot(x)*/
225 t1 = t2 = s = 1.;
226 for(k=0, j=k-n; j < m; k++, j++, s = -s) {
227 /* k == n+j , s = (-1)^k */
228 t1 *= M_PI;/* t1 == pi^(k+1) */
229 if(k >= 2)
230 t2 *= k;/* t2 == k! == gamma(k+1) */
231 if(j >= 0) /* by cheat above, tt === d_k(x) */
232 ans[j] = s*(ans[j] + t1/t2 * tt);
233 }
234 if (n == 0 && kode == 2) /* unused from R, but "wrong": xln === 0 :*/
235 ans[0] += xln;
236 return;
237 } /* x <= 0 */
238
239 /* else : x > 0 */
240 *nz = 0;
241 xln = log(x);
242 if(kode == 1 && m == 1) {/* the R case --- for very large x: */
243 double lrg = 1/(2. * DBL_EPSILON);
244 if(n == 0 && x * xln > lrg) {
245 ans[0] = -xln;
246 return;
247 }
248 else if(n >= 1 && x > n * lrg) {
249 ans[0] = exp(-n * xln)/n; /* == x^-n / n == 1/(n * x^n) */
250 return;
251 }
252 }
253 mm = m;
254 nx = imin2(-Rf_i1mach(15), Rf_i1mach(16));/* = 1021 */
255 r1m5 = Rf_d1mach(5);
256 r1m4 = Rf_d1mach(4) * 0.5;
257 wdtol = fmax2(r1m4, 0.5e-18); /* 1.11e-16 */
258
259 /* elim = approximate exponential over and underflow limit */
260 elim = 2.302 * (nx * r1m5 - 3.0);/* = 700.6174... */
261 for(;;) {
262 nn = n + mm - 1;
263 fn = nn;
264 t = (fn + 1) * xln;
265
266 /* overflow and underflow test for small and large x */
267
268 if (fabs(t) > elim) {
269 if (t <= 0.0) {
270 *nz = 0;
271 *ierr = 2;
272 return;
273 }
274 }
275 else {
276 if (x < wdtol) {
277 ans[0] = pow(x, -n-1.0);
278 if (mm != 1) {
279 for(k = 1; k < mm ; k++)
280 ans[k] = ans[k-1] / x;
281 }
282 if (n == 0 && kode == 2)
283 ans[0] += xln;
284 return;
285 }
286
287 /* compute xmin and the number of terms of the series, fln+1 */
288
289 rln = r1m5 * Rf_i1mach(14);
290 rln = fmin2(rln, 18.06);
291 fln = fmax2(rln, 3.0) - 3.0;
292 yint = 3.50 + 0.40 * fln;
293 slope = 0.21 + fln * (0.0006038 * fln + 0.008677);
294 xm = yint + slope * fn;
295 mx = (int)xm + 1;
296 xmin = mx;
297 if (n != 0) {
298 xm = -2.302 * rln - fmin2(0.0, xln);
299 arg = xm / n;
300 arg = fmin2(0.0, arg);
301 eps = exp(arg);
302 xm = 1.0 - eps;
303 if (fabs(arg) < 1.0e-3)
304 xm = -arg;
305 fln = x * xm / eps;
306 xm = xmin - x;
307 if (xm > 7.0 && fln < 15.0)
308 break;
309 }
310 xdmy = x;
311 xdmln = xln;
312 xinc = 0.0;
313 if (x < xmin) {
314 nx = (int)x;
315 xinc = xmin - nx;
316 xdmy = x + xinc;
317 xdmln = log(xdmy);
318 }
319
320 /* generate w(n+mm-1, x) by the asymptotic expansion */
321
322 t = fn * xdmln;
323 t1 = xdmln + xdmln;
324 t2 = t + xdmln;
325 tk = fmax2(fabs(t), fmax2(fabs(t1), fabs(t2)));
326 if (tk <= elim) /* for all but large x */
327 goto L10;
328 }
329 nz++; /* underflow */
330 mm--;
331 ans[mm] = 0.;
332 if (mm == 0)
333 return;
334 } /* end{for()} */
335 nn = (int)fln + 1;
336 np = n + 1;
337 t1 = (n + 1) * xln;
338 t = exp(-t1);
339 s = t;
340 den = x;
341 for(i=1; i <= nn; i++) {
342 den += 1.;
343 trm[i] = pow(den, (double)-np);
344 s += trm[i];
345 }
346 ans[0] = s;
347 if (n == 0 && kode == 2)
348 ans[0] = s + xln;
349
350 if (mm != 1) { /* generate higher derivatives, j > n */
351
352 tol = wdtol / 5.0;
353 for(j = 1; j < mm; j++) {
354 t /= x;
355 s = t;
356 tols = t * tol;
357 den = x;
358 for(i=1; i <= nn; i++) {
359 den += 1.;
360 trm[i] /= den;
361 s += trm[i];
362 if (trm[i] < tols)
363 break;
364 }
365 ans[j] = s;
366 }
367 }
368 return;
369
370 L10:
371 tss = exp(-t);
372 tt = 0.5 / xdmy;
373 t1 = tt;
374 tst = wdtol * tt;
375 if (nn != 0)
376 t1 = tt + 1.0 / fn;
377 rxsq = 1.0 / (xdmy * xdmy);
378 ta = 0.5 * rxsq;
379 t = (fn + 1) * ta;
380 s = t * bvalues[2];
381 if (fabs(s) >= tst) {
382 tk = 2.0;
383 for(k = 4; k <= 22; k++) {
384 t = t * ((tk + fn + 1)/(tk + 1.0))*((tk + fn)/(tk + 2.0)) * rxsq;
385 trm[k] = t * bvalues[k-1];
386 if (fabs(trm[k]) < tst)
387 break;
388 s += trm[k];
389 tk += 2.;
390 }
391 }
392 s = (s + t1) * tss;
393 if (xinc != 0.0) {
394
395 /* backward recur from xdmy to x */
396
397 nx = (int)xinc;
398 np = nn + 1;
399 if (nx > n_max) {
400 *nz = 0;
401 *ierr = 3;
402 return;
403 }
404 else {
405 if (nn==0)
406 goto L20;
407 xm = xinc - 1.0;
408 fx = x + xm;
409
410 /* this loop should not be changed. fx is accurate when x is small */
411 for(i = 1; i <= nx; i++) {
412 trmr[i] = pow(fx, (double)-np);
413 s += trmr[i];
414 xm -= 1.;
415 fx = x + xm;
416 }
417 }
418 }
419 ans[mm-1] = s;
420 if (fn == 0)
421 goto L30;
422
423 /* generate lower derivatives, j < n+mm-1 */
424
425 for(j = 2; j <= mm; j++) {
426 fn--;
427 tss *= xdmy;
428 t1 = tt;
429 if (fn!=0)
430 t1 = tt + 1.0 / fn;
431 t = (fn + 1) * ta;
432 s = t * bvalues[2];
433 if (fabs(s) >= tst) {
434 tk = 4 + fn;
435 for(k=4; k <= 22; k++) {
436 trm[k] = trm[k] * (fn + 1) / tk;
437 if (fabs(trm[k]) < tst)
438 break;
439 s += trm[k];
440 tk += 2.;
441 }
442 }
443 s = (s + t1) * tss;
444 if (xinc != 0.0) {
445 if (fn == 0)
446 goto L20;
447 xm = xinc - 1.0;
448 fx = x + xm;
449 for(i=1 ; i<=nx ; i++) {
450 trmr[i] = trmr[i] * fx;
451 s += trmr[i];
452 xm -= 1.;
453 fx = x + xm;
454 }
455 }
456 ans[mm - j] = s;
457 if (fn == 0)
458 goto L30;
459 }
460 return;
461
462 L20:
463 for(i = 1; i <= nx; i++)
464 s += 1. / (x + (nx - i)); /* avoid disastrous cancellation, PR#13714 */
465
466 L30:
467 if (kode != 2) /* always */
468 ans[0] = s - xdmln;
469 else if (xdmy != x) {
470 xq = xdmy / x;
471 ans[0] = s - log(xq);
472 }
473 return;
474 } /* dpsifn() */
475
476 #ifdef MATHLIB_STANDALONE
477 # define ML_TREAT_psigam(_IERR_) \
478 if(_IERR_ != 0) { \
479 errno = EDOM; \
480 return ML_NAN; \
481 }
482 #else
483 # define ML_TREAT_psigam(_IERR_) \
484 if(_IERR_ != 0) \
485 return ML_NAN
486 #endif
487
psigamma(double x,double deriv)488 double psigamma(double x, double deriv)
489 {
490 /* n-th derivative of psi(x); e.g., psigamma(x,0) == digamma(x) */
491 double ans;
492 int nz, ierr, k, n;
493
494 if(ISNAN(x))
495 return x;
496 deriv = floor(deriv + 0.5);
497 n = (int)deriv;
498 if(n > n_max) {
499 MATHLIB_WARNING2(_("deriv = %d > %d (= n_max)\n"), n, n_max);
500 return ML_NAN;
501 }
502 dpsifn(x, n, 1, 1, &ans, &nz, &ierr);
503 ML_TREAT_psigam(ierr);
504 /* Now, ans == A := (-1)^(n+1) / gamma(n+1) * psi(n, x) */
505 ans = -ans; /* = (-1)^(0+1) * gamma(0+1) * A */
506 for(k = 1; k <= n; k++)
507 ans *= (-k);/* = (-1)^(k+1) * gamma(k+1) * A */
508 return ans;/* = psi(n, x) */
509 }
510
digamma(double x)511 double digamma(double x)
512 {
513 double ans;
514 int nz, ierr;
515 if(ISNAN(x)) return x;
516 dpsifn(x, 0, 1, 1, &ans, &nz, &ierr);
517 ML_TREAT_psigam(ierr);
518 return -ans;
519 }
520
trigamma(double x)521 double trigamma(double x)
522 {
523 double ans;
524 int nz, ierr;
525 if(ISNAN(x)) return x;
526 dpsifn(x, 1, 1, 1, &ans, &nz, &ierr);
527 ML_TREAT_psigam(ierr);
528 return ans;
529 }
530
tetragamma(double x)531 double tetragamma(double x)
532 {
533 double ans;
534 int nz, ierr;
535 if(ISNAN(x)) return x;
536 dpsifn(x, 2, 1, 1, &ans, &nz, &ierr);
537 ML_TREAT_psigam(ierr);
538 return -2.0 * ans;
539 }
540
pentagamma(double x)541 double pentagamma(double x)
542 {
543 double ans;
544 int nz, ierr;
545 if(ISNAN(x)) return x;
546 dpsifn(x, 3, 1, 1, &ans, &nz, &ierr);
547 ML_TREAT_psigam(ierr);
548 return 6.0 * ans;
549 }
550