1 /* specfunc/psi.c
2 *
3 * Copyright (C) 2007 Brian Gough
4 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2005, 2006 Gerard Jungman
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3 of the License, or (at
9 * your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19 */
20
21 /* Author: G. Jungman */
22
23 #include "gsl__config.h"
24 #include "gsl_math.h"
25 #include "gsl_errno.h"
26 #include "gsl_sf_exp.h"
27 #include "gsl_sf_gamma.h"
28 #include "gsl_sf_zeta.h"
29 #include "gsl_sf_psi.h"
30 #include "gsl_complex_math.h"
31
32 #include <stdio.h>
33
34 #include "gsl_specfunc__error.h"
35
36 #include "gsl_specfunc__chebyshev.h"
37 #include "gsl_specfunc__cheb_eval.c"
38
39 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
40
41
42 /* Chebyshev fit for f(y) = Re(Psi(1+Iy)) + M_EULER - y^2/(1+y^2) - y^2/(2(4+y^2))
43 * 1 < y < 10
44 * ==>
45 * y(x) = (9x + 11)/2, -1 < x < 1
46 * x(y) = (2y - 11)/9
47 *
48 * g(x) := f(y(x))
49 */
50 static double r1py_data[] = {
51 1.59888328244976954803168395603,
52 0.67905625353213463845115658455,
53 -0.068485802980122530009506482524,
54 -0.005788184183095866792008831182,
55 0.008511258167108615980419855648,
56 -0.004042656134699693434334556409,
57 0.001352328406159402601778462956,
58 -0.000311646563930660566674525382,
59 0.000018507563785249135437219139,
60 0.000028348705427529850296492146,
61 -0.000019487536014574535567541960,
62 8.0709788710834469408621587335e-06,
63 -2.2983564321340518037060346561e-06,
64 3.0506629599604749843855962658e-07,
65 1.3042238632418364610774284846e-07,
66 -1.2308657181048950589464690208e-07,
67 5.7710855710682427240667414345e-08,
68 -1.8275559342450963966092636354e-08,
69 3.1020471300626589420759518930e-09,
70 6.8989327480593812470039430640e-10,
71 -8.7182290258923059852334818997e-10,
72 4.4069147710243611798213548777e-10,
73 -1.4727311099198535963467200277e-10,
74 2.7589682523262644748825844248e-11,
75 4.1871826756975856411554363568e-12,
76 -6.5673460487260087541400767340e-12,
77 3.4487900886723214020103638000e-12,
78 -1.1807251417448690607973794078e-12,
79 2.3798314343969589258709315574e-13,
80 2.1663630410818831824259465821e-15
81 };
82 static cheb_series r1py_cs = {
83 r1py_data,
84 29,
85 -1,1,
86 18
87 };
88
89
90 /* Chebyshev fits from SLATEC code for psi(x)
91
92 Series for PSI on the interval 0. to 1.00000D+00
93 with weighted error 2.03E-17
94 log weighted error 16.69
95 significant figures required 16.39
96 decimal places required 17.37
97
98 Series for APSI on the interval 0. to 2.50000D-01
99 with weighted error 5.54E-17
100 log weighted error 16.26
101 significant figures required 14.42
102 decimal places required 16.86
103
104 */
105
106 static double psics_data[23] = {
107 -.038057080835217922,
108 .491415393029387130,
109 -.056815747821244730,
110 .008357821225914313,
111 -.001333232857994342,
112 .000220313287069308,
113 -.000037040238178456,
114 .000006283793654854,
115 -.000001071263908506,
116 .000000183128394654,
117 -.000000031353509361,
118 .000000005372808776,
119 -.000000000921168141,
120 .000000000157981265,
121 -.000000000027098646,
122 .000000000004648722,
123 -.000000000000797527,
124 .000000000000136827,
125 -.000000000000023475,
126 .000000000000004027,
127 -.000000000000000691,
128 .000000000000000118,
129 -.000000000000000020
130 };
131 static cheb_series psi_cs = {
132 psics_data,
133 22,
134 -1, 1,
135 17
136 };
137
138 static double apsics_data[16] = {
139 -.0204749044678185,
140 -.0101801271534859,
141 .0000559718725387,
142 -.0000012917176570,
143 .0000000572858606,
144 -.0000000038213539,
145 .0000000003397434,
146 -.0000000000374838,
147 .0000000000048990,
148 -.0000000000007344,
149 .0000000000001233,
150 -.0000000000000228,
151 .0000000000000045,
152 -.0000000000000009,
153 .0000000000000002,
154 -.0000000000000000
155 };
156 static cheb_series apsi_cs = {
157 apsics_data,
158 15,
159 -1, 1,
160 9
161 };
162
163 #define PSI_TABLE_NMAX 100
164 static double psi_table[PSI_TABLE_NMAX+1] = {
165 0.0, /* Infinity */ /* psi(0) */
166 -M_EULER, /* psi(1) */
167 0.42278433509846713939348790992, /* ... */
168 0.92278433509846713939348790992,
169 1.25611766843180047272682124325,
170 1.50611766843180047272682124325,
171 1.70611766843180047272682124325,
172 1.87278433509846713939348790992,
173 2.01564147795560999653634505277,
174 2.14064147795560999653634505277,
175 2.25175258906672110764745616389,
176 2.35175258906672110764745616389,
177 2.44266167997581201673836525479,
178 2.52599501330914535007169858813,
179 2.60291809023222227314862166505,
180 2.67434666166079370172005023648,
181 2.74101332832746036838671690315,
182 2.80351332832746036838671690315,
183 2.86233685773922507426906984432,
184 2.91789241329478062982462539988,
185 2.97052399224214905087725697883,
186 3.02052399224214905087725697883,
187 3.06814303986119666992487602645,
188 3.11359758531574212447033057190,
189 3.15707584618530734186163491973,
190 3.1987425128519740085283015864,
191 3.2387425128519740085283015864,
192 3.2772040513135124700667631249,
193 3.3142410883505495071038001619,
194 3.3499553740648352213895144476,
195 3.3844381326855248765619282407,
196 3.4177714660188582098952615740,
197 3.4500295305349872421533260902,
198 3.4812795305349872421533260902,
199 3.5115825608380175451836291205,
200 3.5409943255438998981248055911,
201 3.5695657541153284695533770196,
202 3.5973435318931062473311547974,
203 3.6243705589201332743581818244,
204 3.6506863483938174848844976139,
205 3.6763273740348431259101386396,
206 3.7013273740348431259101386396,
207 3.7257176179372821503003825420,
208 3.7495271417468059598241920658,
209 3.7727829557002943319172153216,
210 3.7955102284275670591899425943,
211 3.8177324506497892814121648166,
212 3.8394715810845718901078169905,
213 3.8607481768292527411716467777,
214 3.8815815101625860745049801110,
215 3.9019896734278921969539597029,
216 3.9219896734278921969539597029,
217 3.9415975165651470989147440166,
218 3.9608282857959163296839747858,
219 3.9796962103242182164764276160,
220 3.9982147288427367349949461345,
221 4.0163965470245549168131279527,
222 4.0342536898816977739559850956,
223 4.0517975495308205809735289552,
224 4.0690389288411654085597358518,
225 4.0859880813835382899156680552,
226 4.1026547480502049565823347218,
227 4.1190481906731557762544658694,
228 4.1351772229312202923834981274,
229 4.1510502388042361653993711433,
230 4.1666752388042361653993711433,
231 4.1820598541888515500147557587,
232 4.1972113693403667015299072739,
233 4.2121367424746950597388624977,
234 4.2268426248276362362094507330,
235 4.2413353784508246420065521823,
236 4.2556210927365389277208378966,
237 4.2697055997787924488475984600,
238 4.2835944886676813377364873489,
239 4.2972931188046676391063503626,
240 4.3108066323181811526198638761,
241 4.3241399656515144859531972094,
242 4.3372978603883565912163551041,
243 4.3502848733753695782293421171,
244 4.3631053861958823987421626300,
245 4.3757636140439836645649474401,
246 4.3882636140439836645649474401,
247 4.4006092930563293435772931191,
248 4.4128044150075488557724150703,
249 4.4248526077786331931218126607,
250 4.4367573696833950978837174226,
251 4.4485220755657480390601880108,
252 4.4601499825424922251066996387,
253 4.4716442354160554434975042364,
254 4.4830078717796918071338678728,
255 4.4942438268358715824147667492,
256 4.5053549379469826935258778603,
257 4.5163439489359936825368668713,
258 4.5272135141533849868846929582,
259 4.5379662023254279976373811303,
260 4.5486045001977684231692960239,
261 4.5591308159872421073798223397,
262 4.5695474826539087740464890064,
263 4.5798567610044242379640147796,
264 4.5900608426370772991885045755,
265 4.6001618527380874001986055856
266 };
267
268
269 #define PSI_1_TABLE_NMAX 100
270 static double psi_1_table[PSI_1_TABLE_NMAX+1] = {
271 0.0, /* Infinity */ /* psi(1,0) */
272 M_PI*M_PI/6.0, /* psi(1,1) */
273 0.644934066848226436472415, /* ... */
274 0.394934066848226436472415,
275 0.2838229557371153253613041,
276 0.2213229557371153253613041,
277 0.1813229557371153253613041,
278 0.1535451779593375475835263,
279 0.1331370146940314251345467,
280 0.1175120146940314251345467,
281 0.1051663356816857461222010,
282 0.0951663356816857461222010,
283 0.0869018728717683907503002,
284 0.0799574284273239463058557,
285 0.0740402686640103368384001,
286 0.0689382278476838062261552,
287 0.0644937834032393617817108,
288 0.0605875334032393617817108,
289 0.0571273257907826143768665,
290 0.0540409060376961946237801,
291 0.0512708229352031198315363,
292 0.0487708229352031198315363,
293 0.0465032492390579951149830,
294 0.0444371335365786562720078,
295 0.0425467743683366902984728,
296 0.0408106632572255791873617,
297 0.0392106632572255791873617,
298 0.0377313733163971768204978,
299 0.0363596312039143235969038,
300 0.0350841209998326909438426,
301 0.0338950603577399442137594,
302 0.0327839492466288331026483,
303 0.0317433665203020901265817,
304 0.03076680402030209012658168,
305 0.02984853037475571730748159,
306 0.02898347847164153045627052,
307 0.02816715194102928555831133,
308 0.02739554700275768062003973,
309 0.02666508681283803124093089,
310 0.02597256603721476254286995,
311 0.02531510384129102815759710,
312 0.02469010384129102815759710,
313 0.02409521984367056414807896,
314 0.02352832641963428296894063,
315 0.02298749353699501850166102,
316 0.02247096461137518379091722,
317 0.02197713745088135663042339,
318 0.02150454765882086513703965,
319 0.02105185413233829383780923,
320 0.02061782635456051606003145,
321 0.02020133322669712580597065,
322 0.01980133322669712580597065,
323 0.01941686571420193164987683,
324 0.01904704322899483105816086,
325 0.01869104465298913508094477,
326 0.01834810912486842177504628,
327 0.01801753061247172756017024,
328 0.01769865306145131939690494,
329 0.01739086605006319997554452,
330 0.01709360088954001329302371,
331 0.01680632711763538818529605,
332 0.01652854933985761040751827,
333 0.01625980437882562975715546,
334 0.01599965869724394401313881,
335 0.01574770606433893015574400,
336 0.01550356543933893015574400,
337 0.01526687904880638577704578,
338 0.01503731063741979257227076,
339 0.01481454387422086185273411,
340 0.01459828089844231513993134,
341 0.01438824099085987447620523,
342 0.01418415935820681325171544,
343 0.01398578601958352422176106,
344 0.01379288478501562298719316,
345 0.01360523231738567365335942,
346 0.01342261726990576130858221,
347 0.01324483949212798353080444,
348 0.01307170929822216635628920,
349 0.01290304679189732236910755,
350 0.01273868124291638877278934,
351 0.01257845051066194236996928,
352 0.01242220051066194236996928,
353 0.01226978472038606978956995,
354 0.01212106372098095378719041,
355 0.01197590477193174490346273,
356 0.01183418141592267460867815,
357 0.01169577311142440471248438,
358 0.01156056489076458859566448,
359 0.01142844704164317229232189,
360 0.01129931481023821361463594,
361 0.01117306812421372175754719,
362 0.01104961133409026496742374,
363 0.01092885297157366069257770,
364 0.01081070552355853781923177,
365 0.01069508522063334415522437,
366 0.01058191183901270133041676,
367 0.01047110851491297833872701,
368 0.01036260157046853389428257,
369 0.01025632035036012704977199, /* ... */
370 0.01015219706839427948625679, /* psi(1,99) */
371 0.01005016666333357139524567 /* psi(1,100) */
372 };
373
374
375 /* digamma for x both positive and negative; we do both
376 * cases here because of the way we use even/odd parts
377 * of the function
378 */
379 static int
psi_x(const double x,gsl_sf_result * result)380 psi_x(const double x, gsl_sf_result * result)
381 {
382 const double y = fabs(x);
383
384 if(x == 0.0 || x == -1.0 || x == -2.0) {
385 DOMAIN_ERROR(result);
386 }
387 else if(y >= 2.0) {
388 const double t = 8.0/(y*y)-1.0;
389 gsl_sf_result result_c;
390 cheb_eval_e(&apsi_cs, t, &result_c);
391 if(x < 0.0) {
392 const double s = sin(M_PI*x);
393 const double c = cos(M_PI*x);
394 if(fabs(s) < 2.0*GSL_SQRT_DBL_MIN) {
395 DOMAIN_ERROR(result);
396 }
397 else {
398 result->val = log(y) - 0.5/x + result_c.val - M_PI * c/s;
399 result->err = M_PI*fabs(x)*GSL_DBL_EPSILON/(s*s);
400 result->err += result_c.err;
401 result->err += GSL_DBL_EPSILON * fabs(result->val);
402 return GSL_SUCCESS;
403 }
404 }
405 else {
406 result->val = log(y) - 0.5/x + result_c.val;
407 result->err = result_c.err;
408 result->err += GSL_DBL_EPSILON * fabs(result->val);
409 return GSL_SUCCESS;
410 }
411 }
412 else { /* -2 < x < 2 */
413 gsl_sf_result result_c;
414
415 if(x < -1.0) { /* x = -2 + v */
416 const double v = x + 2.0;
417 const double t1 = 1.0/x;
418 const double t2 = 1.0/(x+1.0);
419 const double t3 = 1.0/v;
420 cheb_eval_e(&psi_cs, 2.0*v-1.0, &result_c);
421
422 result->val = -(t1 + t2 + t3) + result_c.val;
423 result->err = GSL_DBL_EPSILON * (fabs(t1) + fabs(x/(t2*t2)) + fabs(x/(t3*t3)));
424 result->err += result_c.err;
425 result->err += GSL_DBL_EPSILON * fabs(result->val);
426 return GSL_SUCCESS;
427 }
428 else if(x < 0.0) { /* x = -1 + v */
429 const double v = x + 1.0;
430 const double t1 = 1.0/x;
431 const double t2 = 1.0/v;
432 cheb_eval_e(&psi_cs, 2.0*v-1.0, &result_c);
433
434 result->val = -(t1 + t2) + result_c.val;
435 result->err = GSL_DBL_EPSILON * (fabs(t1) + fabs(x/(t2*t2)));
436 result->err += result_c.err;
437 result->err += GSL_DBL_EPSILON * fabs(result->val);
438 return GSL_SUCCESS;
439 }
440 else if(x < 1.0) { /* x = v */
441 const double t1 = 1.0/x;
442 cheb_eval_e(&psi_cs, 2.0*x-1.0, &result_c);
443
444 result->val = -t1 + result_c.val;
445 result->err = GSL_DBL_EPSILON * t1;
446 result->err += result_c.err;
447 result->err += GSL_DBL_EPSILON * fabs(result->val);
448 return GSL_SUCCESS;
449 }
450 else { /* x = 1 + v */
451 const double v = x - 1.0;
452 return cheb_eval_e(&psi_cs, 2.0*v-1.0, result);
453 }
454 }
455 }
456
457
458 /* psi(z) for large |z| in the right half-plane; [Abramowitz + Stegun, 6.3.18] */
459 static
460 gsl_complex
psi_complex_asymp(gsl_complex z)461 psi_complex_asymp(gsl_complex z)
462 {
463 /* coefficients in the asymptotic expansion for large z;
464 * let w = z^(-2) and write the expression in the form
465 *
466 * ln(z) - 1/(2z) - 1/12 w (1 + c1 w + c2 w + c3 w + ... )
467 */
468 static const double c1 = -0.1;
469 static const double c2 = 1.0/21.0;
470 static const double c3 = -0.05;
471
472 gsl_complex zi = gsl_complex_inverse(z);
473 gsl_complex w = gsl_complex_mul(zi, zi);
474 gsl_complex cs;
475
476 /* Horner method evaluation of term in parentheses */
477 gsl_complex sum;
478 sum = gsl_complex_mul_real(w, c3/c2);
479 sum = gsl_complex_add_real(sum, 1.0);
480 sum = gsl_complex_mul_real(sum, c2/c1);
481 sum = gsl_complex_mul(sum, w);
482 sum = gsl_complex_add_real(sum, 1.0);
483 sum = gsl_complex_mul_real(sum, c1);
484 sum = gsl_complex_mul(sum, w);
485 sum = gsl_complex_add_real(sum, 1.0);
486
487 /* correction added to log(z) */
488 cs = gsl_complex_mul(sum, w);
489 cs = gsl_complex_mul_real(cs, -1.0/12.0);
490 cs = gsl_complex_add(cs, gsl_complex_mul_real(zi, -0.5));
491
492 return gsl_complex_add(gsl_complex_log(z), cs);
493 }
494
495
496
497 /* psi(z) for complex z in the right half-plane */
498 static int
psi_complex_rhp(gsl_complex z,gsl_sf_result * result_re,gsl_sf_result * result_im)499 psi_complex_rhp(
500 gsl_complex z,
501 gsl_sf_result * result_re,
502 gsl_sf_result * result_im
503 )
504 {
505 int n_recurse = 0;
506 int i;
507 gsl_complex a;
508
509 if(GSL_REAL(z) == 0.0 && GSL_IMAG(z) == 0.0)
510 {
511 result_re->val = 0.0;
512 result_im->val = 0.0;
513 result_re->err = 0.0;
514 result_im->err = 0.0;
515 return GSL_EDOM;
516 }
517
518 /* compute the number of recurrences to apply */
519 if(GSL_REAL(z) < 20.0 && fabs(GSL_IMAG(z)) < 20.0)
520 {
521 const double sp = sqrt(20.0 + GSL_IMAG(z));
522 const double sn = sqrt(20.0 - GSL_IMAG(z));
523 const double rhs = sp*sn - GSL_REAL(z);
524 if(rhs > 0.0) n_recurse = ceil(rhs);
525 }
526
527 /* compute asymptotic at the large value z + n_recurse */
528 a = psi_complex_asymp(gsl_complex_add_real(z, n_recurse));
529
530 result_re->err = 2.0 * GSL_DBL_EPSILON * fabs(GSL_REAL(a));
531 result_im->err = 2.0 * GSL_DBL_EPSILON * fabs(GSL_IMAG(a));
532
533 /* descend recursively, if necessary */
534 for(i = n_recurse; i >= 1; --i)
535 {
536 gsl_complex zn = gsl_complex_add_real(z, i - 1.0);
537 gsl_complex zn_inverse = gsl_complex_inverse(zn);
538 a = gsl_complex_sub(a, zn_inverse);
539
540 /* accumulate the error, to catch cancellations */
541 result_re->err += 2.0 * GSL_DBL_EPSILON * fabs(GSL_REAL(zn_inverse));
542 result_im->err += 2.0 * GSL_DBL_EPSILON * fabs(GSL_IMAG(zn_inverse));
543 }
544
545 result_re->val = GSL_REAL(a);
546 result_im->val = GSL_IMAG(a);
547
548 result_re->err += 2.0 * GSL_DBL_EPSILON * fabs(result_re->val);
549 result_im->err += 2.0 * GSL_DBL_EPSILON * fabs(result_im->val);
550
551 return GSL_SUCCESS;
552 }
553
554
555
556 /* generic polygamma; assumes n >= 0 and x > 0
557 */
558 static int
psi_n_xg0(const int n,const double x,gsl_sf_result * result)559 psi_n_xg0(const int n, const double x, gsl_sf_result * result)
560 {
561 if(n == 0) {
562 return gsl_sf_psi_e(x, result);
563 }
564 else {
565 /* Abramowitz + Stegun 6.4.10 */
566 gsl_sf_result ln_nf;
567 gsl_sf_result hzeta;
568 int stat_hz = gsl_sf_hzeta_e(n+1.0, x, &hzeta);
569 int stat_nf = gsl_sf_lnfact_e((unsigned int) n, &ln_nf);
570 int stat_e = gsl_sf_exp_mult_err_e(ln_nf.val, ln_nf.err,
571 hzeta.val, hzeta.err,
572 result);
573 if(GSL_IS_EVEN(n)) result->val = -result->val;
574 return GSL_ERROR_SELECT_3(stat_e, stat_nf, stat_hz);
575 }
576 }
577
578
579
580 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
581
gsl_sf_psi_int_e(const int n,gsl_sf_result * result)582 int gsl_sf_psi_int_e(const int n, gsl_sf_result * result)
583 {
584 /* CHECK_POINTER(result) */
585
586 if(n <= 0) {
587 DOMAIN_ERROR(result);
588 }
589 else if(n <= PSI_TABLE_NMAX) {
590 result->val = psi_table[n];
591 result->err = GSL_DBL_EPSILON * fabs(result->val);
592 return GSL_SUCCESS;
593 }
594 else {
595 /* Abramowitz+Stegun 6.3.18 */
596 const double c2 = -1.0/12.0;
597 const double c3 = 1.0/120.0;
598 const double c4 = -1.0/252.0;
599 const double c5 = 1.0/240.0;
600 const double ni2 = (1.0/n)*(1.0/n);
601 const double ser = ni2 * (c2 + ni2 * (c3 + ni2 * (c4 + ni2*c5)));
602 result->val = log(n) - 0.5/n + ser;
603 result->err = GSL_DBL_EPSILON * (fabs(log(n)) + fabs(0.5/n) + fabs(ser));
604 result->err += GSL_DBL_EPSILON * fabs(result->val);
605 return GSL_SUCCESS;
606 }
607 }
608
609
gsl_sf_psi_e(const double x,gsl_sf_result * result)610 int gsl_sf_psi_e(const double x, gsl_sf_result * result)
611 {
612 /* CHECK_POINTER(result) */
613 return psi_x(x, result);
614 }
615
616
617 int
gsl_sf_psi_1piy_e(const double y,gsl_sf_result * result)618 gsl_sf_psi_1piy_e(const double y, gsl_sf_result * result)
619 {
620 const double ay = fabs(y);
621
622 /* CHECK_POINTER(result) */
623
624 if(ay > 1000.0) {
625 /* [Abramowitz+Stegun, 6.3.19] */
626 const double yi2 = 1.0/(ay*ay);
627 const double lny = log(ay);
628 const double sum = yi2 * (1.0/12.0 + 1.0/120.0 * yi2 + 1.0/252.0 * yi2*yi2);
629 result->val = lny + sum;
630 result->err = 2.0 * GSL_DBL_EPSILON * (fabs(lny) + fabs(sum));
631 return GSL_SUCCESS;
632 }
633 else if(ay > 10.0) {
634 /* [Abramowitz+Stegun, 6.3.19] */
635 const double yi2 = 1.0/(ay*ay);
636 const double lny = log(ay);
637 const double sum = yi2 * (1.0/12.0 +
638 yi2 * (1.0/120.0 +
639 yi2 * (1.0/252.0 +
640 yi2 * (1.0/240.0 +
641 yi2 * (1.0/132.0 + 691.0/32760.0 * yi2)))));
642 result->val = lny + sum;
643 result->err = 2.0 * GSL_DBL_EPSILON * (fabs(lny) + fabs(sum));
644 return GSL_SUCCESS;
645 }
646 else if(ay > 1.0){
647 const double y2 = ay*ay;
648 const double x = (2.0*ay - 11.0)/9.0;
649 const double v = y2*(1.0/(1.0+y2) + 0.5/(4.0+y2));
650 gsl_sf_result result_c;
651 cheb_eval_e(&r1py_cs, x, &result_c);
652 result->val = result_c.val - M_EULER + v;
653 result->err = result_c.err;
654 result->err += 2.0 * GSL_DBL_EPSILON * (fabs(v) + M_EULER + fabs(result_c.val));
655 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
656 result->err *= 5.0; /* FIXME: losing a digit somewhere... maybe at x=... ? */
657 return GSL_SUCCESS;
658 }
659 else {
660 /* [Abramowitz+Stegun, 6.3.17]
661 *
662 * -M_EULER + y^2 Sum[1/n 1/(n^2 + y^2), {n,1,M}]
663 * + Sum[1/n^3, {n,M+1,Infinity}]
664 * - y^2 Sum[1/n^5, {n,M+1,Infinity}]
665 * + y^4 Sum[1/n^7, {n,M+1,Infinity}]
666 * - y^6 Sum[1/n^9, {n,M+1,Infinity}]
667 * + O(y^8)
668 *
669 * We take M=50 for at least 15 digit precision.
670 */
671 const int M = 50;
672 const double y2 = y*y;
673 const double c0 = 0.00019603999466879846570;
674 const double c2 = 3.8426659205114376860e-08;
675 const double c4 = 1.0041592839497643554e-11;
676 const double c6 = 2.9516743763500191289e-15;
677 const double p = c0 + y2 *(-c2 + y2*(c4 - y2*c6));
678 double sum = 0.0;
679 double v;
680
681 int n;
682 for(n=1; n<=M; n++) {
683 sum += 1.0/(n * (n*n + y*y));
684 }
685
686 v = y2 * (sum + p);
687 result->val = -M_EULER + v;
688 result->err = GSL_DBL_EPSILON * (M_EULER + fabs(v));
689 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
690 return GSL_SUCCESS;
691 }
692 }
693
694
gsl_sf_psi_1_int_e(const int n,gsl_sf_result * result)695 int gsl_sf_psi_1_int_e(const int n, gsl_sf_result * result)
696 {
697 /* CHECK_POINTER(result) */
698 if(n <= 0) {
699 DOMAIN_ERROR(result);
700 }
701 else if(n <= PSI_1_TABLE_NMAX) {
702 result->val = psi_1_table[n];
703 result->err = GSL_DBL_EPSILON * result->val;
704 return GSL_SUCCESS;
705 }
706 else {
707 /* Abramowitz+Stegun 6.4.12
708 * double-precision for n > 100
709 */
710 const double c0 = -1.0/30.0;
711 const double c1 = 1.0/42.0;
712 const double c2 = -1.0/30.0;
713 const double ni2 = (1.0/n)*(1.0/n);
714 const double ser = ni2*ni2 * (c0 + ni2*(c1 + c2*ni2));
715 result->val = (1.0 + 0.5/n + 1.0/(6.0*n*n) + ser) / n;
716 result->err = GSL_DBL_EPSILON * result->val;
717 return GSL_SUCCESS;
718 }
719 }
720
721
gsl_sf_psi_1_e(const double x,gsl_sf_result * result)722 int gsl_sf_psi_1_e(const double x, gsl_sf_result * result)
723 {
724 /* CHECK_POINTER(result) */
725
726 if(x == 0.0 || x == -1.0 || x == -2.0) {
727 DOMAIN_ERROR(result);
728 }
729 else if(x > 0.0)
730 {
731 return psi_n_xg0(1, x, result);
732 }
733 else if(x > -5.0)
734 {
735 /* Abramowitz + Stegun 6.4.6 */
736 int M = -floor(x);
737 double fx = x + M;
738 double sum = 0.0;
739 int m;
740
741 if(fx == 0.0)
742 DOMAIN_ERROR(result);
743
744 for(m = 0; m < M; ++m)
745 sum += 1.0/((x+m)*(x+m));
746
747 {
748 int stat_psi = psi_n_xg0(1, fx, result);
749 result->val += sum;
750 result->err += M * GSL_DBL_EPSILON * sum;
751 return stat_psi;
752 }
753 }
754 else
755 {
756 /* Abramowitz + Stegun 6.4.7 */
757 const double sin_px = sin(M_PI * x);
758 const double d = M_PI*M_PI/(sin_px*sin_px);
759 gsl_sf_result r;
760 int stat_psi = psi_n_xg0(1, 1.0-x, &r);
761 result->val = d - r.val;
762 result->err = r.err + 2.0*GSL_DBL_EPSILON*d;
763 return stat_psi;
764 }
765 }
766
767
gsl_sf_psi_n_e(const int n,const double x,gsl_sf_result * result)768 int gsl_sf_psi_n_e(const int n, const double x, gsl_sf_result * result)
769 {
770 /* CHECK_POINTER(result) */
771
772 if(n == 0)
773 {
774 return gsl_sf_psi_e(x, result);
775 }
776 else if(n == 1)
777 {
778 return gsl_sf_psi_1_e(x, result);
779 }
780 else if(n < 0 || x <= 0.0) {
781 DOMAIN_ERROR(result);
782 }
783 else {
784 gsl_sf_result ln_nf;
785 gsl_sf_result hzeta;
786 int stat_hz = gsl_sf_hzeta_e(n+1.0, x, &hzeta);
787 int stat_nf = gsl_sf_lnfact_e((unsigned int) n, &ln_nf);
788 int stat_e = gsl_sf_exp_mult_err_e(ln_nf.val, ln_nf.err,
789 hzeta.val, hzeta.err,
790 result);
791 if(GSL_IS_EVEN(n)) result->val = -result->val;
792 return GSL_ERROR_SELECT_3(stat_e, stat_nf, stat_hz);
793 }
794 }
795
796
797 int
gsl_sf_complex_psi_e(const double x,const double y,gsl_sf_result * result_re,gsl_sf_result * result_im)798 gsl_sf_complex_psi_e(
799 const double x,
800 const double y,
801 gsl_sf_result * result_re,
802 gsl_sf_result * result_im
803 )
804 {
805 if(x >= 0.0)
806 {
807 gsl_complex z = gsl_complex_rect(x, y);
808 return psi_complex_rhp(z, result_re, result_im);
809 }
810 else
811 {
812 /* reflection formula [Abramowitz+Stegun, 6.3.7] */
813 gsl_complex z = gsl_complex_rect(x, y);
814 gsl_complex omz = gsl_complex_rect(1.0 - x, -y);
815 gsl_complex zpi = gsl_complex_mul_real(z, M_PI);
816 gsl_complex cotzpi = gsl_complex_cot(zpi);
817 int ret_val = psi_complex_rhp(omz, result_re, result_im);
818
819 if(GSL_IS_REAL(GSL_REAL(cotzpi)) && GSL_IS_REAL(GSL_IMAG(cotzpi)))
820 {
821 result_re->val -= M_PI * GSL_REAL(cotzpi);
822 result_im->val -= M_PI * GSL_IMAG(cotzpi);
823 return ret_val;
824 }
825 else
826 {
827 GSL_ERROR("singularity", GSL_EDOM);
828 }
829 }
830 }
831
832
833
834 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
835
836 #include "gsl_specfunc__eval.h"
837
gsl_sf_psi_int(const int n)838 double gsl_sf_psi_int(const int n)
839 {
840 EVAL_RESULT(gsl_sf_psi_int_e(n, &result));
841 }
842
gsl_sf_psi(const double x)843 double gsl_sf_psi(const double x)
844 {
845 EVAL_RESULT(gsl_sf_psi_e(x, &result));
846 }
847
gsl_sf_psi_1piy(const double x)848 double gsl_sf_psi_1piy(const double x)
849 {
850 EVAL_RESULT(gsl_sf_psi_1piy_e(x, &result));
851 }
852
gsl_sf_psi_1_int(const int n)853 double gsl_sf_psi_1_int(const int n)
854 {
855 EVAL_RESULT(gsl_sf_psi_1_int_e(n, &result));
856 }
857
gsl_sf_psi_1(const double x)858 double gsl_sf_psi_1(const double x)
859 {
860 EVAL_RESULT(gsl_sf_psi_1_e(x, &result));
861 }
862
gsl_sf_psi_n(const int n,const double x)863 double gsl_sf_psi_n(const int n, const double x)
864 {
865 EVAL_RESULT(gsl_sf_psi_n_e(n, x, &result));
866 }
867