1 #include <gnumeric-config.h>
2 #include <sf-gamma.h>
3 #include <sf-trig.h>
4 #include <mathfunc.h>
5
6 #define ML_ERR_return_NAN { return gnm_nan; }
7 #define ML_UNDERFLOW (GNM_EPSILON * GNM_EPSILON)
8 #define ML_ERROR(cause) do { } while(0)
9
10 static int qgammaf (gnm_float x, GnmQuad *mant, int *exp2);
11 static void pochhammer_small_n (gnm_float x, gnm_float n, GnmQuad *res);
12
13
14 /* Compute gnm_log(gamma(a+1)) accurately also for small a (0 < a < 0.5). */
lgamma1p(gnm_float a)15 gnm_float lgamma1p (gnm_float a)
16 {
17 const gnm_float eulers_const = GNM_const(0.5772156649015328606065120900824024);
18
19 /* coeffs[i] holds (zeta(i+2)-1)/(i+2) , i = 1:N, N = 40 : */
20 const int N = 40;
21 static const gnm_float coeffs[40] = {
22 GNM_const(0.3224670334241132182362075833230126e-0),
23 GNM_const(0.6735230105319809513324605383715000e-1),
24 GNM_const(0.2058080842778454787900092413529198e-1),
25 GNM_const(0.7385551028673985266273097291406834e-2),
26 GNM_const(0.2890510330741523285752988298486755e-2),
27 GNM_const(0.1192753911703260977113935692828109e-2),
28 GNM_const(0.5096695247430424223356548135815582e-3),
29 GNM_const(0.2231547584535793797614188036013401e-3),
30 GNM_const(0.9945751278180853371459589003190170e-4),
31 GNM_const(0.4492623673813314170020750240635786e-4),
32 GNM_const(0.2050721277567069155316650397830591e-4),
33 GNM_const(0.9439488275268395903987425104415055e-5),
34 GNM_const(0.4374866789907487804181793223952411e-5),
35 GNM_const(0.2039215753801366236781900709670839e-5),
36 GNM_const(0.9551412130407419832857179772951265e-6),
37 GNM_const(0.4492469198764566043294290331193655e-6),
38 GNM_const(0.2120718480555466586923135901077628e-6),
39 GNM_const(0.1004322482396809960872083050053344e-6),
40 GNM_const(0.4769810169363980565760193417246730e-7),
41 GNM_const(0.2271109460894316491031998116062124e-7),
42 GNM_const(0.1083865921489695409107491757968159e-7),
43 GNM_const(0.5183475041970046655121248647057669e-8),
44 GNM_const(0.2483674543802478317185008663991718e-8),
45 GNM_const(0.1192140140586091207442548202774640e-8),
46 GNM_const(0.5731367241678862013330194857961011e-9),
47 GNM_const(0.2759522885124233145178149692816341e-9),
48 GNM_const(0.1330476437424448948149715720858008e-9),
49 GNM_const(0.6422964563838100022082448087644648e-10),
50 GNM_const(0.3104424774732227276239215783404066e-10),
51 GNM_const(0.1502138408075414217093301048780668e-10),
52 GNM_const(0.7275974480239079662504549924814047e-11),
53 GNM_const(0.3527742476575915083615072228655483e-11),
54 GNM_const(0.1711991790559617908601084114443031e-11),
55 GNM_const(0.8315385841420284819798357793954418e-12),
56 GNM_const(0.4042200525289440065536008957032895e-12),
57 GNM_const(0.1966475631096616490411045679010286e-12),
58 GNM_const(0.9573630387838555763782200936508615e-13),
59 GNM_const(0.4664076026428374224576492565974577e-13),
60 GNM_const(0.2273736960065972320633279596737272e-13),
61 GNM_const(0.1109139947083452201658320007192334e-13)
62 };
63
64 const gnm_float c = GNM_const(0.2273736845824652515226821577978691e-12);/* zeta(N+2)-1 */
65 gnm_float lgam;
66 int i;
67
68 if (gnm_abs (a) >= 0.5)
69 return gnm_lgamma (a + 1);
70
71 /* Abramowitz & Stegun 6.1.33,
72 * also http://functions.wolfram.com/06.11.06.0008.01 */
73 lgam = c * gnm_logcf (-a / 2, N + 2, 1);
74 for (i = N - 1; i >= 0; i--)
75 lgam = coeffs[i] - a * lgam;
76
77 return (a * lgam - eulers_const) * a - log1pmx (a);
78 } /* lgamma1p */
79
80 /* ------------------------------------------------------------------------ */
81
82 /* Imported src/nmath/stirlerr.c from R. */
83 /*
84 * AUTHOR
85 * Catherine Loader, catherine@research.bell-labs.com.
86 * October 23, 2000.
87 *
88 * Merge in to R:
89 * Copyright (C) 2000, The R Core Development Team
90 *
91 * This program is free software; you can redistribute it and/or modify
92 * it under the terms of the GNU General Public License as published by
93 * the Free Software Foundation; either version 2 of the License, or
94 * (at your option) any later version.
95 *
96 * This program is distributed in the hope that it will be useful,
97 * but WITHOUT ANY WARRANTY; without even the implied warranty of
98 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
99 * GNU General Public License for more details.
100 *
101 * You should have received a copy of the GNU General Public License
102 * along with this program; if not, write to the Free Software
103 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
104 * USA.
105 *
106 *
107 * DESCRIPTION
108 *
109 * Computes the log of the error term in Stirling's formula.
110 * For n > 15, uses the series 1/12n - 1/360n^3 + ...
111 * For n <=15, integers or half-integers, uses stored values.
112 * For other n < 15, uses lgamma directly (don't use this to
113 * write lgamma!)
114 *
115 * Merge in to R:
116 * Copyright (C) 2000, The R Core Development Team
117 * R has lgammafn, and lgamma is not part of ISO C
118 */
119
120
121 /* stirlerr(n) = gnm_log(n!) - gnm_log( gnm_sqrt(2*pi*n)*(n/e)^n )
122 * = gnm_log Gamma(n+1) - 1/2 * [gnm_log(2*pi) + gnm_log(n)] - n*[gnm_log(n) - 1]
123 * = gnm_log Gamma(n+1) - (n + 1/2) * gnm_log(n) + n - gnm_log(2*pi)/2
124 *
125 * see also lgammacor() in ./lgammacor.c which computes almost the same!
126 */
127
stirlerr(gnm_float n)128 gnm_float stirlerr(gnm_float n)
129 {
130
131 #define S0 GNM_const(0.083333333333333333333) /* 1/12 */
132 #define S1 GNM_const(0.00277777777777777777778) /* 1/360 */
133 #define S2 GNM_const(0.00079365079365079365079365) /* 1/1260 */
134 #define S3 GNM_const(0.000595238095238095238095238) /* 1/1680 */
135 #define S4 GNM_const(0.0008417508417508417508417508)/* 1/1188 */
136
137 /*
138 error for 0, 0.5, 1.0, 1.5, ..., 14.5, 15.0.
139 */
140 static const gnm_float sferr_halves[31] = {
141 0.0, /* n=0 - wrong, place holder only */
142 GNM_const(0.1534264097200273452913848), /* 0.5 */
143 GNM_const(0.0810614667953272582196702), /* 1.0 */
144 GNM_const(0.0548141210519176538961390), /* 1.5 */
145 GNM_const(0.0413406959554092940938221), /* 2.0 */
146 GNM_const(0.03316287351993628748511048), /* 2.5 */
147 GNM_const(0.02767792568499833914878929), /* 3.0 */
148 GNM_const(0.02374616365629749597132920), /* 3.5 */
149 GNM_const(0.02079067210376509311152277), /* 4.0 */
150 GNM_const(0.01848845053267318523077934), /* 4.5 */
151 GNM_const(0.01664469118982119216319487), /* 5.0 */
152 GNM_const(0.01513497322191737887351255), /* 5.5 */
153 GNM_const(0.01387612882307074799874573), /* 6.0 */
154 GNM_const(0.01281046524292022692424986), /* 6.5 */
155 GNM_const(0.01189670994589177009505572), /* 7.0 */
156 GNM_const(0.01110455975820691732662991), /* 7.5 */
157 GNM_const(0.010411265261972096497478567), /* 8.0 */
158 GNM_const(0.009799416126158803298389475), /* 8.5 */
159 GNM_const(0.009255462182712732917728637), /* 9.0 */
160 GNM_const(0.008768700134139385462952823), /* 9.5 */
161 GNM_const(0.008330563433362871256469318), /* 10.0 */
162 GNM_const(0.007934114564314020547248100), /* 10.5 */
163 GNM_const(0.007573675487951840794972024), /* 11.0 */
164 GNM_const(0.007244554301320383179543912), /* 11.5 */
165 GNM_const(0.006942840107209529865664152), /* 12.0 */
166 GNM_const(0.006665247032707682442354394), /* 12.5 */
167 GNM_const(0.006408994188004207068439631), /* 13.0 */
168 GNM_const(0.006171712263039457647532867), /* 13.5 */
169 GNM_const(0.005951370112758847735624416), /* 14.0 */
170 GNM_const(0.005746216513010115682023589), /* 14.5 */
171 GNM_const(0.005554733551962801371038690) /* 15.0 */
172 };
173 gnm_float nn;
174
175 if (n <= 15.0) {
176 nn = n + n;
177 if (nn == (int)nn) return(sferr_halves[(int)nn]);
178 return(lgamma1p (n ) - (n + 0.5)*gnm_log(n) + n - M_LN_SQRT_2PI);
179 }
180
181 nn = n*n;
182 if (n>500) return((S0-S1/nn)/n);
183 if (n> 80) return((S0-(S1-S2/nn)/nn)/n);
184 if (n> 35) return((S0-(S1-(S2-S3/nn)/nn)/nn)/n);
185 /* 15 < n <= 35 : */
186 return((S0-(S1-(S2-(S3-S4/nn)/nn)/nn)/nn)/n);
187 }
188 /* Cleaning up done by tools/import-R: */
189 #undef S0
190 #undef S1
191 #undef S2
192 #undef S3
193 #undef S4
194
195 /* ------------------------------------------------------------------------ */
196 /* Imported src/nmath/chebyshev.c from R. */
197 /*
198 * Mathlib : A C Library of Special Functions
199 * Copyright (C) 1998 Ross Ihaka
200 *
201 * This program is free software; you can redistribute it and/or modify
202 * it under the terms of the GNU General Public License as published by
203 * the Free Software Foundation; either version 2 of the License, or
204 * (at your option) any later version.
205 *
206 * This program is distributed in the hope that it will be useful,
207 * but WITHOUT ANY WARRANTY; without even the implied warranty of
208 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
209 * GNU General Public License for more details.
210 *
211 * You should have received a copy of the GNU General Public License
212 * along with this program; if not, write to the Free Software
213 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
214 * 02110-1301 USA.
215 *
216 * SYNOPSIS
217 *
218 * int chebyshev_init(double *dos, int nos, double eta)
219 * double chebyshev_eval(double x, double *a, int n)
220 *
221 * DESCRIPTION
222 *
223 * "chebyshev_init" determines the number of terms for the
224 * double precision orthogonal series "dos" needed to insure
225 * the error is no larger than "eta". Ordinarily eta will be
226 * chosen to be one-tenth machine precision.
227 *
228 * "chebyshev_eval" evaluates the n-term Chebyshev series
229 * "a" at "x".
230 *
231 * NOTES
232 *
233 * These routines are translations into C of Fortran routines
234 * by W. Fullerton of Los Alamos Scientific Laboratory.
235 *
236 * Based on the Fortran routine dcsevl by W. Fullerton.
237 * Adapted from R. Broucke, Algorithm 446, CACM., 16, 254 (1973).
238 */
239
240
241 /* NaNs propagated correctly */
242
243
244 /* Definition of function chebyshev_init removed. */
245
246
chebyshev_eval(gnm_float x,const gnm_float * a,const int n)247 static gnm_float chebyshev_eval(gnm_float x, const gnm_float *a, const int n)
248 {
249 gnm_float b0, b1, b2, twox;
250 int i;
251
252 if (n < 1 || n > 1000) ML_ERR_return_NAN;
253
254 if (x < -1.1 || x > 1.1) ML_ERR_return_NAN;
255
256 twox = x * 2;
257 b2 = b1 = 0;
258 b0 = 0;
259 for (i = 1; i <= n; i++) {
260 b2 = b1;
261 b1 = b0;
262 b0 = twox * b1 - b2 + a[n - i];
263 }
264 return (b0 - b2) * 0.5;
265 }
266
267 /* ------------------------------------------------------------------------ */
268 /* Imported src/nmath/lgammacor.c from R. */
269 /*
270 * Mathlib : A C Library of Special Functions
271 * Copyright (C) 1998 Ross Ihaka
272 * Copyright (C) 2000-2001 The R Development Core Team
273 *
274 * This program is free software; you can redistribute it and/or modify
275 * it under the terms of the GNU General Public License as published by
276 * the Free Software Foundation; either version 2 of the License, or
277 * (at your option) any later version.
278 *
279 * This program is distributed in the hope that it will be useful,
280 * but WITHOUT ANY WARRANTY; without even the implied warranty of
281 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
282 * GNU General Public License for more details.
283 *
284 * You should have received a copy of the GNU General Public License
285 * along with this program; if not, write to the Free Software
286 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
287 * USA.
288 *
289 * SYNOPSIS
290 *
291 * #include <Rmath.h>
292 * double lgammacor(double x);
293 *
294 * DESCRIPTION
295 *
296 * Compute the log gamma correction factor for x >= 10 so that
297 *
298 * log(gamma(x)) = .5*log(2*pi) + (x-.5)*log(x) -x + lgammacor(x)
299 *
300 * [ lgammacor(x) is called Del(x) in other contexts (e.g. dcdflib)]
301 *
302 * NOTES
303 *
304 * This routine is a translation into C of a Fortran subroutine
305 * written by W. Fullerton of Los Alamos Scientific Laboratory.
306 *
307 * SEE ALSO
308 *
309 * Loader(1999)'s stirlerr() {in ./stirlerr.c} is *very* similar in spirit,
310 * is faster and cleaner, but is only defined "fast" for half integers.
311 */
312
313
lgammacor(gnm_float x)314 static gnm_float lgammacor(gnm_float x)
315 {
316 static const gnm_float algmcs[15] = {
317 GNM_const(+.1666389480451863247205729650822e+0),
318 GNM_const(-.1384948176067563840732986059135e-4),
319 GNM_const(+.9810825646924729426157171547487e-8),
320 GNM_const(-.1809129475572494194263306266719e-10),
321 GNM_const(+.6221098041892605227126015543416e-13),
322 GNM_const(-.3399615005417721944303330599666e-15),
323 GNM_const(+.2683181998482698748957538846666e-17),
324 GNM_const(-.2868042435334643284144622399999e-19),
325 GNM_const(+.3962837061046434803679306666666e-21),
326 GNM_const(-.6831888753985766870111999999999e-23),
327 GNM_const(+.1429227355942498147573333333333e-24),
328 GNM_const(-.3547598158101070547199999999999e-26),
329 GNM_const(+.1025680058010470912000000000000e-27),
330 GNM_const(-.3401102254316748799999999999999e-29),
331 GNM_const(+.1276642195630062933333333333333e-30)
332 };
333
334 gnm_float tmp;
335
336 #ifdef NOMORE_FOR_THREADS
337 static int nalgm = 0;
338 static gnm_float xbig = 0, xmax = 0;
339
340 /* Initialize machine dependent constants, the first time gamma() is called.
341 FIXME for threads ! */
342 if (nalgm == 0) {
343 /* For IEEE gnm_float precision : nalgm = 5 */
344 nalgm = chebyshev_init(algmcs, 15, GNM_EPSILON/2);/*was d1mach(3)*/
345 xbig = 1 / gnm_sqrt(GNM_EPSILON/2); /* ~ 94906265.6 for IEEE gnm_float */
346 xmax = gnm_exp(fmin2(gnm_log(GNM_MAX / 12), -gnm_log(12 * GNM_MIN)));
347 /* = GNM_MAX / 48 ~= 3.745e306 for IEEE gnm_float */
348 }
349 #else
350 /* For IEEE gnm_float precision GNM_EPSILON = 2^-52 = GNM_const(2.220446049250313e-16) :
351 * xbig = 2 ^ 26.5
352 * xmax = GNM_MAX / 48 = 2^1020 / 3 */
353 # define nalgm 5
354 # define xbig GNM_const(94906265.62425156)
355 # define xmax GNM_const(3.745194030963158e306)
356 #endif
357
358 if (x < 10)
359 ML_ERR_return_NAN
360 else if (x >= xmax) {
361 ML_ERROR(ME_UNDERFLOW);
362 return ML_UNDERFLOW;
363 }
364 else if (x < xbig) {
365 tmp = 10 / x;
366 return chebyshev_eval(tmp * tmp * 2 - 1, algmcs, nalgm) / x;
367 }
368 else return 1 / (x * 12);
369 }
370 /* Cleaning up done by tools/import-R: */
371 #undef nalgm
372 #undef xbig
373 #undef xmax
374
375 /* ------------------------------------------------------------------------ */
376 /* Imported src/nmath/lbeta.c from R. */
377 /*
378 * Mathlib : A C Library of Special Functions
379 * Copyright (C) 1998 Ross Ihaka
380 * Copyright (C) 2000 The R Development Core Team
381 * Copyright (C) 2003 The R Foundation
382 *
383 * This program is free software; you can redistribute it and/or modify
384 * it under the terms of the GNU General Public License as published by
385 * the Free Software Foundation; either version 2 of the License, or
386 * (at your option) any later version.
387 *
388 * This program is distributed in the hope that it will be useful,
389 * but WITHOUT ANY WARRANTY; without even the implied warranty of
390 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
391 * GNU General Public License for more details.
392 *
393 * You should have received a copy of the GNU General Public License
394 * along with this program; if not, write to the Free Software
395 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
396 * USA.
397 *
398 * SYNOPSIS
399 *
400 * #include <Rmath.h>
401 * double lbeta(double a, double b);
402 *
403 * DESCRIPTION
404 *
405 * This function returns the value of the log beta function.
406 *
407 * NOTES
408 *
409 * This routine is a translation into C of a Fortran subroutine
410 * by W. Fullerton of Los Alamos Scientific Laboratory.
411 */
412
413
gnm_lbeta(gnm_float a,gnm_float b)414 gnm_float gnm_lbeta(gnm_float a, gnm_float b)
415 {
416 gnm_float corr, p, q;
417
418 p = q = a;
419 if(b < p) p = b;/* := min(a,b) */
420 if(b > q) q = b;/* := max(a,b) */
421
422 #ifdef IEEE_754
423 if(gnm_isnan(a) || gnm_isnan(b))
424 return a + b;
425 #endif
426
427 /* both arguments must be >= 0 */
428
429 if (p < 0)
430 ML_ERR_return_NAN
431 else if (p == 0) {
432 return gnm_pinf;
433 }
434 else if (!gnm_finite(q)) {
435 return gnm_ninf;
436 }
437
438 if (p >= 10) {
439 /* p and q are big. */
440 corr = lgammacor(p) + lgammacor(q) - lgammacor(p + q);
441 return gnm_log(q) * -0.5 + M_LN_SQRT_2PI + corr
442 + (p - 0.5) * gnm_log(p / (p + q)) + q * gnm_log1p(-p / (p + q));
443 }
444 else if (q >= 10) {
445 /* p is small, but q is big. */
446 corr = lgammacor(q) - lgammacor(p + q);
447 return gnm_lgamma(p) + corr + p - p * gnm_log(p + q)
448 + (q - 0.5) * gnm_log1p(-p / (p + q));
449 }
450 else
451 /* p and q are small: p <= q < 10. */
452 return gnm_lgamma (p) + gnm_lgamma (q) - gnm_lgamma (p + q);
453 }
454
455 /* ------------------------------------------------------------------------ */
456
457 gnm_float
gnm_gammax(gnm_float x,int * exp2)458 gnm_gammax (gnm_float x, int *exp2)
459 {
460 GnmQuad r;
461 (void) qgammaf (x, &r, exp2);
462 return gnm_quad_value (&r);
463 }
464
465 /**
466 * gnm_gamma:
467 * @x: a number
468 *
469 * Returns: the Gamma function evaluated at @x for positive or
470 * non-integer @x.
471 */
472 gnm_float
gnm_gamma(gnm_float x)473 gnm_gamma (gnm_float x)
474 {
475 int e;
476 gnm_float r = gnm_gammax (x, &e);
477 return gnm_ldexp (r, e);
478 }
479
480 /* ------------------------------------------------------------------------- */
481
482 gnm_float
gnm_factx(gnm_float x,int * exp2)483 gnm_factx (gnm_float x, int *exp2)
484 {
485 GnmQuad r;
486 (void)qfactf (x, &r, exp2);
487 return gnm_quad_value (&r);
488 }
489
490 /**
491 * gnm_fact:
492 * @x: number
493 *
494 * Returns: the factorial of @x, which must not be a negative integer.
495 */
496 gnm_float
gnm_fact(gnm_float x)497 gnm_fact (gnm_float x)
498 {
499 int e;
500 gnm_float r = gnm_factx (x, &e);
501 return gnm_ldexp (r, e);
502 }
503
504 /* ------------------------------------------------------------------------- */
505
506 /* 0: ok, 1: overflow, 2: nan */
507 static int
qbetaf(gnm_float a,gnm_float b,GnmQuad * mant,int * exp2)508 qbetaf (gnm_float a, gnm_float b, GnmQuad *mant, int *exp2)
509 {
510 GnmQuad ma, mb, mab;
511 int ea, eb, eab;
512 gnm_float ab = a + b;
513
514 if (gnm_isnan (ab) ||
515 (a <= 0 && a == gnm_floor (a)) ||
516 (b <= 0 && b == gnm_floor (b)))
517 return 2;
518
519 if (ab <= 0 && ab == gnm_floor (ab)) {
520 gnm_quad_init (mant, 0);
521 *exp2 = 0;
522 return 0;
523 }
524
525 if (b > a) {
526 gnm_float s = a;
527 a = b;
528 b = s;
529 }
530
531 if (a > 1 && gnm_abs (b) < 1) {
532 void *state;
533 if (qgammaf (b, &mb, &eb))
534 return 1;
535 state = gnm_quad_start ();
536 pochhammer_small_n (a, b, &ma);
537 gnm_quad_div (mant, &mb, &ma);
538 gnm_quad_end (state);
539 *exp2 = eb;
540 return 0;
541 }
542
543 if (!qgammaf (a, &ma, &ea) &&
544 !qgammaf (b, &mb, &eb) &&
545 !qgammaf (ab, &mab, &eab)) {
546 void *state = gnm_quad_start ();
547 gnm_quad_mul (&ma, &ma, &mb);
548 gnm_quad_div (mant, &ma, &mab);
549 gnm_quad_end (state);
550 *exp2 = ea + eb - eab;
551 return 0;
552 } else
553 return 1;
554 }
555
556 /**
557 * gnm_beta:
558 * @a: a number
559 * @b: a number
560 *
561 * Returns: the Beta function evaluated at @a and @b.
562 */
563 gnm_float
gnm_beta(gnm_float a,gnm_float b)564 gnm_beta (gnm_float a, gnm_float b)
565 {
566 GnmQuad r;
567 int e;
568
569 switch (qbetaf (a, b, &r, &e)) {
570 case 0: return gnm_ldexp (gnm_quad_value (&r), e);
571 case 1: return gnm_pinf;
572 default: return gnm_nan;
573 }
574 }
575
576 /**
577 * gnm_lbeta3:
578 * @a: a number
579 * @b: a number
580 * @sign: (out): the sign
581 *
582 * Returns: the logarithm of the absolute value of the Beta function
583 * evaluated at @a and @b. The sign will be stored in @sign as -1 or
584 * +1. This function is useful because the result of the beta
585 * function can be too large for doubles.
586 */
587 gnm_float
gnm_lbeta3(gnm_float a,gnm_float b,int * sign)588 gnm_lbeta3 (gnm_float a, gnm_float b, int *sign)
589 {
590 int sign_a, sign_b, sign_ab;
591 gnm_float ab = a + b;
592 gnm_float res_a, res_b, res_ab;
593 GnmQuad r;
594 int e;
595
596 switch (qbetaf (a, b, &r, &e)) {
597 case 0: {
598 gnm_float m = gnm_quad_value (&r);
599 *sign = (m >= 0 ? +1 : -1);
600 return gnm_log (gnm_abs (m)) + e * M_LN2gnum;
601 }
602 case 1:
603 /* Overflow */
604 break;
605 default:
606 *sign = 1;
607 return gnm_nan;
608 }
609
610 if (a > 0 && b > 0) {
611 *sign = 1;
612 return gnm_lbeta (a, b);
613 }
614
615 /* This is awful */
616 res_a = gnm_lgamma_r (a, &sign_a);
617 res_b = gnm_lgamma_r (b, &sign_b);
618 res_ab = gnm_lgamma_r (ab, &sign_ab);
619
620 *sign = sign_a * sign_b * sign_ab;
621 return res_a + res_b - res_ab;
622 }
623
624 /* ------------------------------------------------------------------------- */
625 /*
626 * This computes the E(x) such that
627 *
628 * Gamma(x) = sqrt(2Pi) * x^(x-1/2) * exp(-x) * E(x)
629 *
630 * x should be >0 and the result is, roughly, 1+1/(12x).
631 */
632 static void
gamma_error_factor(GnmQuad * res,const GnmQuad * x)633 gamma_error_factor (GnmQuad *res, const GnmQuad *x)
634 {
635 gnm_float num[] = {
636 GNM_const(1.),
637 GNM_const(1.),
638 GNM_const(-139.),
639 GNM_const(-571.),
640 GNM_const(163879.),
641 GNM_const(5246819.),
642 GNM_const(-534703531.),
643 GNM_const(-4483131259.),
644 GNM_const(432261921612371.)
645 };
646 gnm_float den[] = {
647 GNM_const(12.),
648 GNM_const(288.),
649 GNM_const(51840.),
650 GNM_const(2488320.),
651 GNM_const(209018880.),
652 GNM_const(75246796800.),
653 GNM_const(902961561600.),
654 GNM_const(86684309913600.),
655 GNM_const(514904800886784000.)
656 };
657 GnmQuad zn, xpn;
658 int i;
659 gnm_float xval = gnm_quad_value (x);
660 int n;
661
662 g_return_if_fail (xval > 0);
663
664 // We want x >= 20 for the asymptotic expansion
665 n = (xval < 20) ? (int)gnm_floor (21 - xval) : 0;
666 gnm_quad_init (&xpn, n);
667 gnm_quad_add (&xpn, &xpn, x);
668
669 gnm_quad_init (&zn, 1);
670 *res = zn;
671
672 for (i = 0; i < (int)G_N_ELEMENTS (num); i++) {
673 GnmQuad t, c;
674 gnm_quad_mul (&zn, &zn, &xpn);
675 gnm_quad_init (&c, den[i]);
676 gnm_quad_mul (&t, &zn, &c);
677 gnm_quad_init (&c, num[i]);
678 gnm_quad_div (&t, &c, &t);
679 gnm_quad_add (res, res, &t);
680 }
681
682 if (n > 0) {
683 int i;
684 GnmQuad en, xxn, xph;
685
686 // Gamma(x) = sqrt(2Pi) * x^(x-1/2) * exp(-x) * E(x)
687 // Gamma(x+n) = sqrt(2Pi) * (x+n)^(x+n-1/2) * exp(-x-n) * E(x+n)
688
689 // E(x+n) / E(x) =
690 // Gamma(x+n)/Gamma(x) * (x^(x-1/2) * exp(-x)) / ((x+n)^(x+n-1/2) * exp(-x-n)) =
691 // (x*(x+1)*...*(x+n-1)) * exp(n) * (x/(x+n))^(x-1/2) / (x+n)^n =
692 // ((x+1)*...*(x+n-1)) * exp(n) * (x/(x+n))^(x+1/2) / (x+n)^(n-1) =
693 // ((x+1)/(x+n)*...*(x+n-1)/(x+n)) * exp(n) * (x/(x+n))^(x+1/2)
694
695 for (i = 1; i < n; i++) {
696 // *= (x+i)/(x+n)
697 GnmQuad xpi;
698 gnm_quad_init (&xpi, i);
699 gnm_quad_add (&xpi, &xpi, x);
700 gnm_quad_div (res, res, &xpi);
701 gnm_quad_mul (res, res, &xpn);
702 }
703
704 // /= exp(n)
705 gnm_quad_init (&en, n);
706 gnm_quad_exp (&en, NULL, &en);
707 gnm_quad_div (res, res, &en);
708
709 // /= (x/(x+n))^(x+1/2)
710 gnm_quad_init (&xph, 0.5);
711 gnm_quad_add (&xph, &xph, x);
712 gnm_quad_div (&xxn, x, &xpn);
713 gnm_quad_pow (&xxn, NULL, &xxn, &xph);
714 gnm_quad_div (res, res, &xxn);
715 }
716 }
717
718 /* ------------------------------------------------------------------------- */
719
720 static void
pochhammer_small_n(gnm_float x,gnm_float n,GnmQuad * res)721 pochhammer_small_n (gnm_float x, gnm_float n, GnmQuad *res)
722 {
723 GnmQuad qx, qn, qr, qs, f1, f2, f3, f4, f5;
724 gnm_float r;
725 gboolean debug = FALSE;
726
727 g_return_if_fail (x >= 1);
728 g_return_if_fail (gnm_abs (n) <= 1);
729
730 /*
731 * G(x) = c * x^(x-1/2) * exp(-x) * E(x)
732 * G(x+n) = c * (x+n)^(x+n-1/2) * exp(-(x+n)) * E(x+n)
733 * = c * (x+n)^(x-1/2) * (x+n)^n * exp(-x) * exp(-n) * E(x+n)
734 *
735 * G(x+n)/G(x)
736 * = (1+n/x)^(x-1/2) * (x+n)^n * exp(-n) * E(x+n)/E(x)
737 * = (1+n/x)^x / sqrt(1+n/x) * (x+n)^n * exp(-n) * E(x+n)/E(x)
738 * = exp(x*log(1+n/x) - n) / sqrt(1+n/x) * (x+n)^n * E(x+n)/E(x)
739 * = exp(x*log1p(n/x) - n) / sqrt(1+n/x) * (x+n)^n * E(x+n)/E(x)
740 * = exp(x*(log1pmx(n/x)+n/x) - n) / sqrt(1+n/x) * (x+n)^n * E(x+n)/E(x)
741 * = exp(x*log1pmx(n/x) + n - n) / sqrt(1+n/x) * (x+n)^n * E(x+n)/E(x)
742 * = exp(x*log1pmx(n/x)) / sqrt(1+n/x) * (x+n)^n * E(x+n)/E(x)
743 */
744
745 gnm_quad_init (&qx, x);
746 gnm_quad_init (&qn, n);
747
748 gnm_quad_div (&qr, &qn, &qx);
749 r = gnm_quad_value (&qr);
750
751 gnm_quad_add (&qs, &qx, &qn);
752
753 /* exp(x*log1pmx(n/x)) */
754 gnm_quad_mul12 (&f1, log1pmx (r), x); /* sub-opt */
755 gnm_quad_exp (&f1, NULL, &f1);
756 if (debug) g_printerr ("f1=%.20g\n", gnm_quad_value (&f1));
757
758 /* sqrt(1+n/x) */
759 gnm_quad_add (&f2, &gnm_quad_one, &qr);
760 gnm_quad_sqrt (&f2, &f2);
761 if (debug) g_printerr ("f2=%.20g\n", gnm_quad_value (&f2));
762
763 /* (x+n)^n */
764 gnm_quad_pow (&f3, NULL, &qs, &qn);
765 if (debug) g_printerr ("f3=%.20g\n", gnm_quad_value (&f3));
766
767 /* E(x+n) */
768 gamma_error_factor (&f4, &qs);
769 if (debug) g_printerr ("f4=%.20g\n", gnm_quad_value (&f4));
770
771 /* E(x) */
772 gamma_error_factor (&f5, &qx);
773 if (debug) g_printerr ("f5=%.20g\n", gnm_quad_value (&f5));
774
775 gnm_quad_div (res, &f1, &f2);
776 gnm_quad_mul (res, res, &f3);
777 gnm_quad_mul (res, res, &f4);
778 gnm_quad_div (res, res, &f5);
779 }
780
781 static gnm_float
pochhammer_naive(gnm_float x,int n)782 pochhammer_naive (gnm_float x, int n)
783 {
784 void *state = gnm_quad_start ();
785 GnmQuad qp, qx;
786 gnm_float r;
787
788 qp = gnm_quad_one;
789 gnm_quad_init (&qx, x);
790 while (n-- > 0) {
791 gnm_quad_mul (&qp, &qp, &qx);
792 gnm_quad_add (&qx, &qx, &gnm_quad_one);
793 }
794 r = gnm_quad_value (&qp);
795 gnm_quad_end (state);
796
797 return r;
798 }
799
800 /**
801 * pochhammer:
802 * @x: a real number
803 * @n: a real number, often an integer
804 *
805 * This function computes Pochhammer's symbol at @x and @n, i.e.,
806 * Gamma(@x+@n)/Gamma(@x). This is well defined unless @x or @x+@n is a
807 * non-negative integer. The ratio has a removable singlularity at @n=0
808 * and the result is 1.
809 *
810 * Returns: Pochhammer's symbol (@x)_@n.
811 */
812 gnm_float
pochhammer(gnm_float x,gnm_float n)813 pochhammer (gnm_float x, gnm_float n)
814 {
815 gnm_float rn, rx, lr;
816 GnmQuad m1, m2;
817 int e1, e2;
818
819 if (gnm_isnan (x) || gnm_isnan (n))
820 return gnm_nan;
821
822 if (n == 0)
823 return 1;
824
825 rx = gnm_floor (x);
826 rn = gnm_floor (n);
827
828 /*
829 * Use naive multiplication when n is a small integer.
830 * We don't want to use this if x is also an integer
831 * (but we might do so below if x is insanely large).
832 */
833 if (n == rn && x != rx && n >= 0 && n < 40)
834 return pochhammer_naive (x, (int)n);
835
836 if (!qfactf (x + n - 1, &m1, &e1) &&
837 !qfactf (x - 1, &m2, &e2)) {
838 void *state = gnm_quad_start ();
839 int de = e1 - e2;
840 GnmQuad qr;
841 gnm_float r;
842
843 gnm_quad_div (&qr, &m1, &m2);
844 r = gnm_quad_value (&qr);
845 gnm_quad_end (state);
846
847 return gnm_ldexp (r, de);
848 }
849
850 if (x == rx && x <= 0) {
851 if (n != rn)
852 return 0;
853 if (x == 0)
854 return (n > 0)
855 ? 0
856 : ((gnm_fmod (-n, 2) == 0 ? +1 : -1) /
857 gnm_fact (-n));
858 if (n > -x)
859 return gnm_nan;
860 }
861
862 /*
863 * We have left the common cases. One of x+n and x is
864 * insanely big, possibly both.
865 */
866
867 if (gnm_abs (x) < 1)
868 return gnm_pinf;
869
870 if (n < 0)
871 return 1 / pochhammer (x + n, -n);
872
873 if (n == rn && n >= 0 && n < 100)
874 return pochhammer_naive (x, (int)n);
875
876 if (gnm_abs (n) < 1) {
877 /* x is big. */
878 void *state = gnm_quad_start ();
879 GnmQuad qr;
880 gnm_float r;
881 pochhammer_small_n (x, n, &qr);
882 r = gnm_quad_value (&qr);
883 gnm_quad_end (state);
884 return r;
885 }
886
887 /* Panic mode. */
888 g_printerr ("x=%.20g n=%.20g\n", x, n);
889 lr = ((x - 0.5) * gnm_log1p (n / x) +
890 n * gnm_log (x + n) -
891 n +
892 (lgammacor (x + n) - lgammacor (x)));
893 return gnm_exp (lr);
894 }
895
896 /* ------------------------------------------------------------------------- */
897
898 static void
rescale_mant_exp(GnmQuad * mant,int * exp2)899 rescale_mant_exp (GnmQuad *mant, int *exp2)
900 {
901 GnmQuad s;
902 int e;
903
904 (void)gnm_frexp (gnm_quad_value (mant), &e);
905 *exp2 += e;
906 gnm_quad_init (&s, gnm_ldexp (1.0, -e));
907 gnm_quad_mul (mant, mant, &s);
908 }
909
910 /* Tabulate up to, but not including, this number. */
911 #define QFACTI_LIMIT 10000
912
913 static gboolean
qfacti(int n,GnmQuad * mant,int * exp2)914 qfacti (int n, GnmQuad *mant, int *exp2)
915 {
916 static GnmQuad mants[QFACTI_LIMIT];
917 static int exp2s[QFACTI_LIMIT];
918 static int init = 0;
919
920 if (n < 0 || n >= QFACTI_LIMIT) {
921 *mant = gnm_quad_zero;
922 *exp2 = 0;
923 return TRUE;
924 }
925
926 if (n >= init) {
927 void *state = gnm_quad_start ();
928
929 if (init == 0) {
930 gnm_quad_init (&mants[0], 0.5);
931 exp2s[0] = 1;
932 init++;
933 }
934
935 while (n >= init) {
936 GnmQuad m;
937
938 gnm_quad_init (&m, init);
939 gnm_quad_mul (&mants[init], &m, &mants[init - 1]);
940 exp2s[init] = exp2s[init - 1];
941 rescale_mant_exp (&mants[init], &exp2s[init]);
942
943 init++;
944 }
945
946 gnm_quad_end (state);
947 }
948
949 *mant = mants[n];
950 *exp2 = exp2s[n];
951 return FALSE;
952 }
953
954 /* 0: ok, 1: overflow, 2: nan */
955 int
qfactf(gnm_float x,GnmQuad * mant,int * exp2)956 qfactf (gnm_float x, GnmQuad *mant, int *exp2)
957 {
958 void *state;
959 gboolean res = 0;
960
961 *exp2 = 0;
962
963 if (gnm_isnan (x) || (x < 0 && x == gnm_floor (x))) {
964 mant->h = mant->l = gnm_nan;
965 return 2;
966 }
967
968 if (x >= G_MAXINT / 2) {
969 mant->h = mant->l = gnm_pinf;
970 return 1;
971 }
972
973 if (x == gnm_floor (x)) {
974 /* 0, 1, 2, ... */
975 if (!qfacti ((int)x, mant, exp2))
976 return 0;
977 }
978
979 state = gnm_quad_start ();
980
981 if (x < -1) {
982 if (qfactf (-x - 1, mant, exp2))
983 res = 1;
984 else {
985 GnmQuad b;
986
987 gnm_quad_init (&b, -x);
988 gnm_quad_sinpi (&b, &b);
989 gnm_quad_mul (&b, &b, mant);
990 gnm_quad_div (mant, &gnm_quad_pi, &b);
991 *exp2 = -*exp2;
992 }
993 } else if (x >= QFACTI_LIMIT - 0.5) {
994 /*
995 * Let y = x + 1 = m * 2^e; c = sqrt(2Pi).
996 *
997 * G(y) = c * y^(y-1/2) * exp(-y) * E(y)
998 * = c * (y/e)^y / sqrt(y) * E(y)
999 */
1000 GnmQuad y, f1, f2, f3, f4;
1001 gnm_float ef2;
1002 gboolean debug = FALSE;
1003
1004 if (debug) g_printerr ("x=%.20g\n", x);
1005
1006 gnm_quad_init (&y, x + 1);
1007 *exp2 = 0;
1008
1009 /* sqrt(2Pi) */
1010 gnm_quad_sqrt (&f1, &gnm_quad_2pi);
1011 if (debug) g_printerr ("f1=%.20g\n", gnm_quad_value (&f1));
1012
1013 /* (y/e)^y */
1014 gnm_quad_div (&f2, &y, &gnm_quad_e);
1015 gnm_quad_pow (&f2, &ef2, &f2, &y);
1016 if (ef2 > G_MAXINT || ef2 < G_MININT) {
1017 res = 1;
1018 f2.h = f2.l = gnm_pinf;
1019 } else
1020 *exp2 += (int)ef2;
1021 if (debug) g_printerr ("f2=%.20g\n", gnm_quad_value (&f2));
1022
1023 /* sqrt(y) */
1024 gnm_quad_sqrt (&f3, &y);
1025 if (debug) g_printerr ("f3=%.20g\n", gnm_quad_value (&f3));
1026
1027 /* E(x) */
1028 gamma_error_factor (&f4, &y);
1029 if (debug) g_printerr ("f4=%.20g\n", gnm_quad_value (&f4));
1030
1031 gnm_quad_mul (mant, &f1, &f2);
1032 gnm_quad_div (mant, mant, &f3);
1033 gnm_quad_mul (mant, mant, &f4);
1034
1035 if (debug) g_printerr ("G(x+1)=%.20g * 2^%d %s\n", gnm_quad_value (mant), *exp2, res ? "overflow" : "");
1036 } else {
1037 GnmQuad s, qx, mFw;
1038 gnm_float w, f;
1039 int eFw;
1040
1041 /*
1042 * w integer, |f|<=0.5, x=w+f.
1043 *
1044 * Do this before we do the stepping below which would kill
1045 * up to 4 bits of accuracy of f.
1046 */
1047 w = gnm_floor (x + 0.5);
1048 f = x - w;
1049 gnm_quad_init (&qx, x);
1050
1051 gnm_quad_init (&s, 1);
1052 while (w < 20) {
1053 gnm_quad_add (&qx, &qx, &gnm_quad_one);
1054 w++;
1055 gnm_quad_mul (&s, &s, &qx);
1056 }
1057
1058 if (qfacti ((int)w, &mFw, &eFw)) {
1059 mant->h = mant->l = gnm_pinf;
1060 res = 1;
1061 } else {
1062 GnmQuad r;
1063
1064 pochhammer_small_n (w + 1, f, &r);
1065 gnm_quad_mul (mant, &mFw, &r);
1066 gnm_quad_div (mant, mant, &s);
1067 *exp2 = eFw;
1068 }
1069 }
1070
1071 if (res == 0)
1072 rescale_mant_exp (mant, exp2);
1073
1074 gnm_quad_end (state);
1075 return res;
1076 }
1077
1078 /* 0: ok, 1: overflow, 2: nan */
1079 static int
qgammaf(gnm_float x,GnmQuad * mant,int * exp2)1080 qgammaf (gnm_float x, GnmQuad *mant, int *exp2)
1081 {
1082 if (x < -1.5 || x > 0.5)
1083 return qfactf (x - 1, mant, exp2);
1084 else if (gnm_isnan (x) || x == gnm_floor (x)) {
1085 *exp2 = 0;
1086 mant->h = mant->l = gnm_nan;
1087 return 2;
1088 } else {
1089 void *state = gnm_quad_start ();
1090 GnmQuad qx;
1091
1092 qfactf (x, mant, exp2);
1093 gnm_quad_init (&qx, x);
1094 gnm_quad_div (mant, mant, &qx);
1095 rescale_mant_exp (mant, exp2);
1096 gnm_quad_end (state);
1097 return 0;
1098 }
1099 }
1100
1101 /* ------------------------------------------------------------------------- */
1102
1103 gnm_float
combin(gnm_float n,gnm_float k)1104 combin (gnm_float n, gnm_float k)
1105 {
1106 GnmQuad m1, m2, m3;
1107 int e1, e2, e3;
1108 gboolean ok;
1109
1110 if (k < 0 || k > n || n != gnm_floor (n) || k != gnm_floor (k))
1111 return gnm_nan;
1112
1113 k = MIN (k, n - k);
1114 if (k == 0)
1115 return 1;
1116 if (k == 1)
1117 return n;
1118
1119 ok = (n < G_MAXINT &&
1120 !qfactf (n, &m1, &e1) &&
1121 !qfactf (k, &m2, &e2) &&
1122 !qfactf (n - k, &m3, &e3));
1123
1124 if (ok) {
1125 void *state = gnm_quad_start ();
1126 gnm_float c;
1127 gnm_quad_mul (&m2, &m2, &m3);
1128 gnm_quad_div (&m1, &m1, &m2);
1129 c = gnm_ldexp (gnm_quad_value (&m1), e1 - e2 - e3);
1130 gnm_quad_end (state);
1131 return c;
1132 }
1133
1134 if (k < 100) {
1135 void *state = gnm_quad_start ();
1136 GnmQuad p, a, b;
1137 gnm_float c;
1138 int i;
1139
1140 gnm_quad_init (&p, 1);
1141 for (i = 0; i < k; i++) {
1142 gnm_quad_init (&a, n - i);
1143 gnm_quad_mul (&p, &p, &a);
1144
1145 gnm_quad_init (&b, i + 1);
1146 gnm_quad_div (&p, &p, &b);
1147 }
1148
1149 c = gnm_quad_value (&p);
1150 gnm_quad_end (state);
1151 return c;
1152 }
1153
1154 return pochhammer (n - k + 1, k) / gnm_fact (k);
1155 }
1156
1157 gnm_float
permut(gnm_float n,gnm_float k)1158 permut (gnm_float n, gnm_float k)
1159 {
1160 if (k < 0 || k > n || n != gnm_floor (n) || k != gnm_floor (k))
1161 return gnm_nan;
1162
1163 return pochhammer (n - k + 1, k);
1164 }
1165
1166 /* ------------------------------------------------------------------------- */
1167
1168 #ifdef GNM_SUPPLIES_LGAMMA
1169 /* Avoid using signgam. It may be missing in system libraries. */
1170 int signgam;
1171
1172 double
lgamma(double x)1173 lgamma (double x)
1174 {
1175 return lgamma_r (x, &signgam);
1176 }
1177 #endif
1178
1179 #ifdef GNM_SUPPLIES_LGAMMA_R
1180 double
lgamma_r(double x,int * signp)1181 lgamma_r (double x, int *signp)
1182 {
1183 *signp = +1;
1184
1185 if (gnm_isnan (x))
1186 return gnm_nan;
1187
1188 if (x > 0) {
1189 gnm_float f = 1;
1190
1191 while (x < 10) {
1192 f *= x;
1193 x++;
1194 }
1195
1196 return (M_LN_SQRT_2PI + (x - 0.5) * gnm_log(x) -
1197 x + lgammacor(x)) - gnm_log (f);
1198 } else {
1199 gnm_float axm2 = gnm_fmod (-x, 2.0);
1200 gnm_float y = gnm_sinpi (axm2) / M_PIgnum;
1201
1202 *signp = axm2 > 1.0 ? +1 : -1;
1203
1204 return y == 0
1205 ? gnm_nan
1206 : - gnm_log (gnm_abs (y)) - lgamma1p (-x);
1207 }
1208 }
1209 #endif
1210
1211 /* ------------------------------------------------------------------------- */
1212
1213
1214 static const gnm_float lanczos_g = GNM_const (808618867.0) / 134217728;
1215
1216 /*
1217 * This Mathematica sniplet computes the Lanczos gamma coefficients:
1218 *
1219 * Dr[k_]:=DiagonalMatrix[Join[{1},Array[-Binomial[2*#-1,#]*#&,k]]]
1220 * c[k_]:= Array[If[#1+#2==2,1/2,If[#1>=#2,(-1)^(#1+#2)*4^(#2-1)*(#1-1)*(#1+#2-3)!/(#1-#2)!/(2*#2-2)!,0]]&,{k+1,k+1}]
1221 * Dc[k_]:=DiagonalMatrix[Array[2*(2*#-3)!!&,k+1]]
1222 * B[k_]:=Array[If[#1==1,1,If[#1<=#2,(-1)^(#2-#1)*Binomial[#1+#2-3,#1+#1-3],0]]&,{k+1,k+1}]
1223 * M[k_]:=(Dr[k].B[k]).(c[k].Dc[k])
1224 * f[g_,k_]:=Array[Sqrt[2]*(E/(2*(#-1+g)+1))^(#-1/2)&,k+1]
1225 * a[g_,k_]:=M[k].f[g,k]*Exp[g]/Sqrt[2*Pi]
1226 *
1227 * The result of a[g,k] will contain both positive and negative constants.
1228 * Most people using the Lanczos series do not understand that a naive
1229 * implementation will suffer significant cancellation errors. The error
1230 * estimates assume the series is computed without loss!
1231 *
1232 * Following Boost we multiply the entire partial fraction by Pochhammer[z+1,k]
1233 * That gives us a polynomium with positive coefficient. For kicks we toss
1234 * the constant factor back in.
1235 *
1236 * b[g_,k_]:=Sum[a[g,k][[i+1]]/If[i==0,1,(z+i)],{i,0,k}]*Pochhammer[z+1,k]
1237 * c13b:=Block[{$MaxExtraPrecision=500},FullSimplify[N[b[808618867/134217728,12]*Sqrt[2*Pi]/Exp[808618867/134217728],300]]]
1238 *
1239 * Finally we recast that in terms of gamma's argument:
1240 *
1241 * N[CoefficientList[c13b /. z->(zp-1),{zp}],50]
1242 *
1243 * Enter complex numbers, exit simplicity. The error bounds for the
1244 * Lanczos approximation are bounds for the absolute value of the result.
1245 * The relative error on one of the coordinates can be much higher.
1246 */
1247 static const gnm_float lanczos_num[] = {
1248 GNM_const(56906521.913471563880907910335591226868592353221448),
1249 GNM_const(103794043.11634454519062710536160702385539539810110),
1250 GNM_const(86363131.288138591455469272889778684223419113014358),
1251 GNM_const(43338889.324676138347737237405905333160850993321475),
1252 GNM_const(14605578.087685068084141699827913592185707234229516),
1253 GNM_const(3481712.1549806459088207101896477455646802362321652),
1254 GNM_const(601859.61716810987866702265336993523025071425740828),
1255 GNM_const(75999.293040145426498753034435989091370919973262979),
1256 GNM_const(6955.9996025153761403563101155151989875259157712039),
1257 GNM_const(449.94455690631681194468586076509884096232715968614),
1258 GNM_const(19.519927882476174828478609662356521362076846583112),
1259 GNM_const(0.50984166556566761881251786448046945099926051133936),
1260 GNM_const(0.0060618423462489065257837539645559368832224636654970)
1261 };
1262
1263 /* CoefficientList[Pochhammer[z,12],z] */
1264 static const guint32 lanczos_denom[G_N_ELEMENTS(lanczos_num)] = {
1265 0, 39916800, 120543840, 150917976, 105258076, 45995730,
1266 13339535, 2637558, 357423, 32670, 1925, 66, 1
1267 };
1268
1269 /**
1270 * gnm_complex_gamma:
1271 * @z: a complex number
1272 * @exp2: (out) (allow-none): Return location for power of 2.
1273 *
1274 * Returns: (transfer full): the Gamma function evaluated at @z.
1275 */
1276 gnm_complex
gnm_complex_gamma(gnm_complex z,int * exp2)1277 gnm_complex_gamma (gnm_complex z, int *exp2)
1278 {
1279 if (exp2)
1280 *exp2 = 0;
1281
1282 if (GNM_CREALP (z)) {
1283 return GNM_CREAL (exp2 ? gnm_gammax (z.re, exp2) : gnm_gamma (z.re));
1284 } else if (z.re < 0) {
1285 /* Gamma(z) = pi / (sin(pi*z) * Gamma(-z+1)) */
1286 gnm_complex b = GNM_CMAKE (M_PIgnum * gnm_fmod (z.re, 2),
1287 M_PIgnum * z.im);
1288 /* Hmm... sin overflows when b.im is large. */
1289 gnm_complex res = GNM_CDIV (GNM_CREAL (M_PIgnum),
1290 GNM_CMUL (gnm_complex_fact (GNM_CNEG (z), exp2),
1291 GNM_CSIN (b)));
1292 if (exp2)
1293 *exp2 = -*exp2;
1294 return res;
1295 } else {
1296 gnm_complex zmh, f, p, q;
1297 int i;
1298
1299 i = G_N_ELEMENTS(lanczos_num) - 1;
1300 p = GNM_CREAL (lanczos_num[i]);
1301 q = GNM_CREAL (lanczos_denom[i]);
1302 while (--i >= 0) {
1303 p = GNM_CMUL (p, z);
1304 p.re += lanczos_num[i];
1305 q = GNM_CMUL (q, z);
1306 q.re += lanczos_denom[i];
1307 }
1308
1309 zmh = GNM_CMAKE (z.re - 0.5, z.im);
1310 f = GNM_CPOW (GNM_CADD (zmh, GNM_CREAL (lanczos_g)),
1311 GNM_CSCALE (zmh, 0.5));
1312
1313 return GNM_CMUL4 (f, GNM_CEXP (GNM_CNEG (zmh)), f, GNM_CDIV (p, q));
1314 }
1315 }
1316
1317 /* ------------------------------------------------------------------------- */
1318
1319 /**
1320 * gnm_complex_fact:
1321 * @z: a complex number
1322 * @exp2: (out) (allow-none): Return location for power of 2.
1323 *
1324 * Returns: (transfer full): the factorial function evaluated at @z.
1325 */
1326 gnm_complex
gnm_complex_fact(gnm_complex z,int * exp2)1327 gnm_complex_fact (gnm_complex z, int *exp2)
1328 {
1329 if (exp2)
1330 *exp2 = 0;
1331
1332 if (GNM_CREALP (z)) {
1333 return GNM_CREAL (exp2 ? gnm_factx (z.re, exp2) : gnm_fact (z.re));
1334 } else {
1335 // This formula is valid for all arguments except zero
1336 // which we conveniently handled above.
1337 return GNM_CMUL (gnm_complex_gamma (z, exp2), z);
1338 }
1339 }
1340
1341 /* ------------------------------------------------------------------------- */
1342
1343 // D(a,z) := z^a * exp(-z) / Gamma (a + 1)
1344 static gnm_complex
complex_temme_D(gnm_complex a,gnm_complex z)1345 complex_temme_D (gnm_complex a, gnm_complex z)
1346 {
1347 gnm_complex t;
1348
1349 // The idea here is to control intermediate sizes and to avoid
1350 // accuracy problems caused by exp and pow. For now, do neither.
1351
1352 t = GNM_CDIV (GNM_CPOW (z, a), GNM_CEXP (z));
1353 return GNM_CDIV (t, gnm_complex_fact (a, NULL));
1354 }
1355
1356
1357 typedef void (*GnmComplexContinuedFraction) (gnm_complex *ai, gnm_complex *bi,
1358 size_t i, const gnm_complex *args);
1359
1360 static gboolean
gnm_complex_continued_fraction(gnm_complex * dst,size_t N,GnmComplexContinuedFraction cf,gnm_complex const * args)1361 gnm_complex_continued_fraction (gnm_complex *dst, size_t N,
1362 GnmComplexContinuedFraction cf,
1363 gnm_complex const *args)
1364 {
1365 gnm_complex A0, A1, B0, B1;
1366 size_t i;
1367 const gboolean debug_cf = FALSE;
1368
1369 A0 = B1 = GNM_C1;
1370 A1 = B0 = GNM_C0;
1371
1372 for (i = 1; i < N; i++) {
1373 gnm_complex ai, bi, t1, t2, c1, c2, A2, B2;
1374 gnm_float m;
1375 const gnm_float BIG = GNM_const(18446744073709551616.0);
1376
1377 cf (&ai, &bi, i, args);
1378
1379 /* Update A. */
1380 t1 = GNM_CMUL (bi, A1);
1381 t2 = GNM_CMUL (ai, A0);
1382 A2 = GNM_CADD (t1, t2);
1383 A0 = A1; A1 = A2;
1384
1385 /* Update B. */
1386 t1 = GNM_CMUL (bi, B1);
1387 t2 = GNM_CMUL (ai, B0);
1388 B2 = GNM_CADD (t1, t2);
1389 B0 = B1; B1 = B2;
1390
1391 /* Rescale */
1392 m = gnm_abs (B1.re) + gnm_abs (B1.im);
1393 if (m >= BIG || m <= 1 / BIG) {
1394 int e;
1395 gnm_float s;
1396
1397 if (m == 0)
1398 return FALSE;
1399
1400 (void)gnm_frexp (m, &e);
1401 if (debug_cf)
1402 g_printerr ("rescale by 2^%d\n", -e);
1403
1404 s = gnm_ldexp (1, -e);
1405 A0 = GNM_CSCALE (A0, s);
1406 A1 = GNM_CSCALE (A1, s);
1407 B0 = GNM_CSCALE (B0, s);
1408 B1 = GNM_CSCALE (B1, s);
1409 }
1410
1411 /* Check for convergence */
1412 t1 = GNM_CMUL (A1, B0);
1413 t2 = GNM_CMUL (A0, B1);
1414 c1 = GNM_CSUB (t1, t2);
1415
1416 c2 = GNM_CMUL (B0, B1);
1417
1418 t1 = GNM_CDIV (A1, B1);
1419 if (debug_cf) {
1420 g_printerr (" a : %.20g + %.20g I\n", ai.re, ai.im);
1421 g_printerr (" b : %.20g + %.20g I\n", bi.re, bi.im);
1422 g_printerr (" A : %.20g + %.20g I\n", A1.re, A1.im);
1423 g_printerr (" B : %.20g + %.20g I\n", B1.re, B1.im);
1424 g_printerr ("%3zd : %.20g + %.20g I\n", i, t1.re, t1.im);
1425 }
1426
1427 if (GNM_CABS (c1) <= GNM_CABS (c2) * (GNM_EPSILON /2))
1428 break;
1429 }
1430
1431 if (i == N) {
1432 g_printerr ("continued fraction failed to converge.\n");
1433 // Make the failure obvious.
1434 *dst = GNM_CNAN;
1435 return FALSE;
1436 }
1437
1438 *dst = GNM_CDIV (A1, B1);
1439 return TRUE;
1440 }
1441
1442 static void
igamma_lower_coefs(gnm_complex * ai,gnm_complex * bi,size_t i,gnm_complex const * args)1443 igamma_lower_coefs (gnm_complex *ai, gnm_complex *bi, size_t i,
1444 gnm_complex const *args)
1445 {
1446 gnm_complex const *a = args + 0;
1447 gnm_complex const *z = args + 1;
1448
1449 if (i == 1)
1450 *ai = GNM_C1;
1451 else if (i & 1) {
1452 *ai = GNM_CSCALE (*z, i >> 1);
1453 } else {
1454 gnm_complex f = GNM_CMAKE (-(a->re + ((i >> 1) - 1)), -a->im);
1455 *ai = GNM_CMUL (f, *z);
1456 }
1457
1458 *bi = GNM_CMAKE (a->re + (i - 1), a->im);
1459 }
1460
1461 static gboolean
igamma_lower_cf(gnm_complex * dst,const gnm_complex * a,const gnm_complex * z)1462 igamma_lower_cf (gnm_complex *dst, const gnm_complex *a, const gnm_complex *z)
1463 {
1464 gnm_complex args[2] = { *a, *z };
1465 gnm_complex res;
1466
1467 if (!gnm_complex_continued_fraction (&res, 100, igamma_lower_coefs, args))
1468 return FALSE;
1469
1470 // FIXME: The following should be handled without creating big numbers.
1471 *dst = GNM_CMUL3 (res, GNM_CEXP (GNM_CNEG (*z)), GNM_CPOW (*z, *a));
1472
1473 return TRUE;
1474 }
1475
1476 static gboolean
igamma_upper_asymp(gnm_complex * dst,const gnm_complex * a,const gnm_complex * z)1477 igamma_upper_asymp (gnm_complex *dst, const gnm_complex *a, const gnm_complex *z)
1478 {
1479 gnm_float am = GNM_CABS (*a);
1480 gnm_float zm = GNM_CABS (*z);
1481 gnm_float n0;
1482 gnm_complex s, t;
1483 gboolean debug = FALSE;
1484 size_t i;
1485
1486 if (am >= zm)
1487 return FALSE;
1488
1489 // Things start to diverge here
1490 n0 = a->re + gnm_sqrt (zm * zm - a->im * a->im);
1491 if (debug)
1492 g_printerr ("n0=%g\n", n0);
1493 (void)n0;
1494
1495 // FIXME: Verify this condition for whether we have enough
1496 // precision at term n0
1497 if (2 * zm < GNM_MANT_DIG * M_LN2gnum)
1498 return FALSE;
1499
1500 s = GNM_C0;
1501
1502 t = complex_temme_D (GNM_CSUB (*a, GNM_C1), *z);
1503
1504 for (i = 0; i < 100; i++) {
1505 s = GNM_CADD (s, t);
1506
1507 if (debug) {
1508 g_printerr ("%3zd: t=%.20g + %.20g * I\n", i, t.re, t.im);
1509 g_printerr (" : s=%.20g + %.20g * I\n", s.re, s.im);
1510 }
1511
1512 if (GNM_CABS (t) <= GNM_CABS (s) * GNM_EPSILON) {
1513 if (debug)
1514 g_printerr ("igamma_upper_asymp converged.\n");
1515 *dst = s;
1516 return TRUE;
1517 }
1518
1519 t = GNM_CDIV (t, *z);
1520 t = GNM_CMUL (t, GNM_CSUB (*a, GNM_CREAL (i + 1)));
1521 }
1522
1523 if (debug)
1524 g_printerr ("igamma_upper_asymp failed to converge.\n");
1525
1526 return FALSE;
1527 }
1528
1529 static void
fixup_upper_real(gnm_complex * res,gnm_complex a,gnm_complex z)1530 fixup_upper_real (gnm_complex *res, gnm_complex a, gnm_complex z)
1531 {
1532 // This assumes we have an upper gamma regularized result.
1533 //
1534 // It appears that upper algorithms have trouble with negative real z
1535 // (for example, such z being outside the allowed domain) in some cases.
1536
1537 if (GNM_CREALP (z) && z.re < 0 &&
1538 GNM_CREALP (a) && a.re != gnm_floor (a.re)) {
1539 // Everything in the lower power series is real except z^a
1540 // which is not. So...
1541 // 1. Flip to lower gamma
1542 // 2. Assume modulus is correct
1543 // 3. Use exact angle for lower gamma
1544 // 4. Flip back to upper gamma
1545 gnm_complex lres = GNM_CSUB (GNM_C1, *res);
1546 *res = GNM_CPOLARPI (GNM_CABS (lres), a.re);
1547 *res = GNM_CSUB (GNM_C1, *res);
1548 }
1549 }
1550
1551 /**
1552 * gnm_complex_igamma:
1553 * @a: a complex number
1554 * @z: a complex number
1555 * @lower: determines if upper or lower incomplete gamma is desired.
1556 * @regularized: determines if the result should be normalized by Gamma(@a).
1557 *
1558 * Returns: (transfer full): the incomplete gamma function evaluated at
1559 * @a and @z.
1560 */
1561 gnm_complex
gnm_complex_igamma(gnm_complex a,gnm_complex z,gboolean lower,gboolean regularized)1562 gnm_complex_igamma (gnm_complex a, gnm_complex z,
1563 gboolean lower, gboolean regularized)
1564 {
1565 gnm_complex res, ga;
1566 gboolean have_lower, have_regularized;
1567 gboolean have_ga = FALSE;
1568
1569 if (regularized && GNM_CREALP (a) &&
1570 a.re <= 0 && a.re == gnm_floor (a.re)) {
1571 res = GNM_C0;
1572 have_lower = FALSE;
1573 have_regularized = TRUE;
1574 goto fixup;
1575 }
1576
1577 if (GNM_CREALP (a) && a.re >= 0 &&
1578 GNM_CREALP (z) && z.re >= 0) {
1579 res = GNM_CREAL (pgamma (z.re, a.re, 1, lower, FALSE));
1580 have_lower = lower;
1581 have_regularized = TRUE;
1582 goto fixup;
1583 }
1584
1585 if (igamma_upper_asymp (&res, &a, &z)) {
1586 have_lower = FALSE;
1587 have_regularized = TRUE;
1588 fixup_upper_real (&res, a, z);
1589 goto fixup;
1590 }
1591
1592 if (igamma_lower_cf (&res, &a, &z)) {
1593 have_lower = TRUE;
1594 have_regularized = FALSE;
1595 goto fixup;
1596 }
1597
1598 // Failure of all sub-methods.
1599 return GNM_CNAN;
1600
1601 fixup:
1602 // Fixup to the desired form as needed. This is not ideal.
1603 // 1. Regularizing here is too late due to overflow.
1604 // 2. Upper/lower switch can have cancellation
1605 if (regularized != have_regularized) {
1606 ga = gnm_complex_gamma (a, NULL);
1607 have_ga = TRUE;
1608
1609 if (have_regularized)
1610 res = GNM_CMUL (res, ga);
1611 else
1612 res = GNM_CDIV (res, ga);
1613 have_regularized = TRUE;
1614 }
1615
1616 if (lower != have_lower) {
1617 if (have_regularized) {
1618 res = GNM_CSUB (GNM_C1, res);
1619 } else {
1620 if (!have_ga)
1621 ga = gnm_complex_gamma (a, NULL);
1622 res = GNM_CSUB (ga, res);
1623 }
1624 }
1625
1626 return res;
1627 }
1628
1629 /* ------------------------------------------------------------------------- */
1630
1631 static gnm_float
gnm_digamma_series_1(gnm_float x)1632 gnm_digamma_series_1 (gnm_float x)
1633 {
1634 static const gnm_float ctr = 3414350731.0 / 4294967296.0; // ~ x0-2/3
1635 // Taylor series data for digamma(x)*x around ctr
1636 // (Multiplying by x eliminates the pole at 0 and improves convergence)
1637
1638 // There are more terms here than will be used in practice
1639 static const gnm_float c[] = {
1640 GNM_const(-1.393604931385844667205297), // cst
1641 GNM_const(+0.7838726021041081531302582), // z
1642 GNM_const(+1.820471535319717826256316), // z^2
1643 GNM_const(+0.2300704039473615371242174), // z^3
1644 GNM_const(-0.03648708728785595477443336), // z^4
1645 GNM_const(+0.008663338335810582341288719), // z^5
1646 GNM_const(-0.002436194723850649598022839), // z^6
1647 GNM_const(+0.0007486622557872255311371203), // z^7
1648 GNM_const(-0.0002423133587459245107417307), // z^8
1649 GNM_const(+0.00008100113830883611703726430), // z^9
1650 GNM_const(-0.00002765115168760370451893173), // z^10
1651 GNM_const(+9.572584786684540889574004e-6), // z^11
1652 GNM_const(-3.345885770126257344664911e-6), // z^12
1653 GNM_const(+1.177300128979172845825083e-6), // z^13
1654 GNM_const(-4.161969426343619044066147e-7), // z^14
1655 GNM_const(+1.476236789046367348142744e-7), // z^15
1656 GNM_const(-5.248645227284800471117817e-8), // z^16
1657 GNM_const(+1.869315129102582931045594e-8), // z^17
1658 GNM_const(-6.665853754670506957976488e-9), // z^18
1659 GNM_const(+2.379136739280974154742874e-9), // z^19
1660 GNM_const(-8.497029470698846358950073e-10), // z^20
1661 GNM_const(+3.036142118559307675850845e-10), // z^21
1662 GNM_const(-1.085246878064370064490199e-10), // z^22
1663 GNM_const(+3.880126402094332901095669e-11), // z^23
1664 GNM_const(-1.387536654151506108828032e-11), // z^24
1665 GNM_const(+4.962524563617018345793409e-12), // z^25
1666 GNM_const(-1.775025683975804156201718e-12), // z^26
1667 GNM_const(+6.349488874733389900536889e-13), // z^27
1668 GNM_const(-2.271415182435993612263917e-13), // z^28
1669 GNM_const(+8.125903477897147090860925e-14), // z^29
1670 GNM_const(-2.907097355266920392577544e-14), // z^30
1671 GNM_const(+1.040056414044071726030447e-14), // z^31
1672 GNM_const(-3.721012573246943428604950e-15), // z^32
1673 GNM_const(+1.331283261904345080561599e-15), // z^33
1674 GNM_const(-4.763032612601286523705145e-16), // z^34
1675 GNM_const(+1.704116960031678756548478e-16), // z^35
1676 GNM_const(-6.097015154289962965498327e-17), // z^36
1677 GNM_const(+2.181406718966981191594648e-17), // z^37
1678 GNM_const(-7.804716283314031896188832e-18), // z^38
1679 GNM_const(+2.792404989185037140120149e-18), // z^39
1680 GNM_const(-9.990800520119412094515637e-19) // z^40
1681 };
1682
1683 gnm_float sum, xn, eps, dx;
1684 unsigned ui;
1685
1686 dx = xn = x - ctr;
1687 sum = c[0] + c[1] * xn;
1688 eps = gnm_abs (GNM_EPSILON / 2 * sum);
1689 for (ui = 2; ui < G_N_ELEMENTS (c); ui++) {
1690 gnm_float t;
1691 xn *= dx;
1692 t = c[ui] * xn;
1693 sum += t;
1694 if (gnm_abs (t) < eps)
1695 break;
1696 }
1697
1698 return sum / x / (x + 1);
1699 }
1700
1701 static gnm_float
gnm_digamma_series_2(gnm_float x,gnm_float dx)1702 gnm_digamma_series_2 (gnm_float x, gnm_float dx)
1703 {
1704 // Taylor series data for digamma(x)*x around x0.
1705 // (Multiplying by x eliminates the pole at 0 and improves convergence)
1706
1707 // There are more terms here than will be used in practice
1708 static const gnm_float c[] = {
1709 0,
1710 GNM_const(1.414380859739958132208209), // z
1711 GNM_const(+0.3205153650531439606356288), // z^2
1712 GNM_const(-0.06493160890417499678330267), // z^3
1713 GNM_const(+0.01887583274794994723362426), // z^4
1714 GNM_const(-0.006343606951359283604253287), // z^5
1715 GNM_const(+0.002294851106215796610898052), // z^6
1716 GNM_const(-0.0008656448634441624396007814), // z^7
1717 GNM_const(+0.0003349197451448133179202073), // z^8
1718 GNM_const(-0.0001316774179498895538138516), // z^9
1719 GNM_const(+0.00005231455693269487786690492), // z^10
1720 GNM_const(-0.00002092930898551028581067484), // z^11
1721 GNM_const(+8.412567983061925868991692e-6), // z^12
1722 GNM_const(-3.392327126020536111624551e-6), // z^13
1723 GNM_const(+1.370973972130058130320036e-6), // z^14
1724 GNM_const(-5.549180707134621401005220e-7), // z^15
1725 GNM_const(+2.248510299244865219544966e-7), // z^16
1726 GNM_const(-9.117750735408115351181446e-8), // z^17
1727 GNM_const(+3.699221275229322519704744e-8), // z^18
1728 GNM_const(-1.501394539608077112213162e-8), // z^19
1729 GNM_const(+6.095280485458728954145874e-9), // z^20
1730 GNM_const(-2.474989843290409518138793e-9), // z^21
1731 GNM_const(+1.005102611341640470198718e-9), // z^22
1732 GNM_const(-4.082140595549856261286344e-10), // z^23
1733 GNM_const(+1.658037290401667848672372e-10), // z^24
1734 GNM_const(-6.734743373121082302361713e-11), // z^25
1735 GNM_const(+2.735661184007454408449954e-11), // z^26
1736 GNM_const(-1.111255343693481217856139e-11), // z^27
1737 GNM_const(+4.514116174157725512713376e-12), // z^28
1738 GNM_const(-1.833735949521707719130688e-12), // z^29
1739 GNM_const(+7.449112972569399411235872e-13), // z^30
1740 GNM_const(-3.026041976266472126189062e-13), // z^31
1741 GNM_const(+1.229269770874759794761121e-13), // z^32
1742 GNM_const(-4.993680849210878859449551e-14), // z^33
1743 GNM_const(+2.028594795343119634764731e-14), // z^34
1744 GNM_const(-8.240821397432176895280819e-15), // z^35
1745 GNM_const(+3.347697235347764148196203e-15), // z^36
1746 GNM_const(-1.359947630194034569577449e-15), // z^37
1747 GNM_const(+5.524569438901596063753430e-16), // z^38
1748 GNM_const(-2.244268739259513574290477e-16), // z^39
1749 GNM_const(+9.116988470680108150341624e-17) // z^40
1750 };
1751 gnm_float sum, xn, eps;
1752 unsigned ui;
1753
1754 xn = dx;
1755 sum = c[1] * xn;
1756 eps = gnm_abs (GNM_EPSILON * sum);
1757 for (ui = 2; ui < G_N_ELEMENTS (c); ui++) {
1758 gnm_float t;
1759 xn *= dx;
1760 t = c[ui] * xn;
1761 sum += t;
1762 if (gnm_abs (t) < eps)
1763 break;
1764 }
1765
1766 return sum / x;
1767 }
1768
1769 static gnm_float
gnm_digamma_series_3(gnm_float x)1770 gnm_digamma_series_3 (gnm_float x)
1771 {
1772 static const gnm_float ctr = 9140973792.0 / 4294967296.0; // ~ x0+2/3
1773
1774 // Taylor series data for digamma(x)*x*(x+1) around ctr.
1775 // (Multiplying by x eliminates the pole at 0 and improves convergence;
1776 // multiplying by x+1 removes trouble caused by the pole at -1.)
1777 //
1778 // There are more terms here than will be used in practice
1779 static const gnm_float c[] = {
1780 GNM_const(1.069187202106379964561108), // cst
1781 GNM_const(+1.772667605096075412537626), // z
1782 GNM_const(+0.2272125634616216308385530), // z^2
1783 GNM_const(-0.03340833758699978856544311), // z^3
1784 GNM_const(+0.007175553429733710899335576), // z^4
1785 GNM_const(-0.001806192980500979068857208), // z^5
1786 GNM_const(+0.0004945960000406938148418368), // z^6
1787 GNM_const(-0.0001423916069504330801643716), // z^7
1788 GNM_const(+0.00004231966722000581929615164), // z^8
1789 GNM_const(-0.00001284637919029649494826060), // z^9
1790 GNM_const(+3.956444268156385801727645e-6), // z^10
1791 GNM_const(-1.230919658902018354780620e-6), // z^11
1792 GNM_const(+3.857326410290438339904409e-7), // z^12
1793 GNM_const(-1.215068812516640310282077e-7), // z^13
1794 GNM_const(+3.842015503145882562666222e-8), // z^14
1795 GNM_const(-1.218213493657190927765958e-8), // z^15
1796 GNM_const(+3.870598142893619365308165e-9), // z^16
1797 GNM_const(-1.231667143855504792729306e-9), // z^17
1798 GNM_const(+3.923744199538871225509428e-10), // z^18
1799 GNM_const(-1.251053017217116281525239e-10), // z^19
1800 GNM_const(+3.991408929102272214329984e-11), // z^20
1801 GNM_const(-1.274041466704381992529912e-11), // z^21
1802 GNM_const(+4.068145278333075741751660e-12), // z^22
1803 GNM_const(-1.299351027914282051289942e-12), // z^23
1804 GNM_const(+4.150924791093867196981568e-13), // z^24
1805 GNM_const(-1.326263739748920828936488e-13), // z^25
1806 GNM_const(+4.238042170781100294818004e-14), // z^26
1807 GNM_const(-1.354374285142548716401069e-14), // z^27
1808 GNM_const(+4.328534597139797982202326e-15), // z^28
1809 GNM_const(-1.383454394822643441972787e-15), // z^29
1810 GNM_const(+4.421862837726160406096067e-16), // z^30
1811 GNM_const(-1.413377420676461469423437e-16), // z^31
1812 GNM_const(+4.517732012277313050103063e-17), // z^32
1813 GNM_const(-1.444075584733824688998366e-17), // z^33
1814 GNM_const(+4.615989316796428759128748e-18), // z^34
1815 GNM_const(-1.475515451985640283943034e-18), // z^35
1816 GNM_const(+4.716565090265976036820537e-19), // z^36
1817 GNM_const(-1.507683748197167538907172e-19), // z^37
1818 GNM_const(+4.819438643661667878185773e-20), // z^38
1819 GNM_const(-1.540579118701962073182260e-20), // z^39
1820 GNM_const(+4.924618369274725064707054e-21) // z^40
1821 };
1822
1823 gnm_float sum, xn, eps, dx;
1824 unsigned ui;
1825
1826 dx = xn = x - ctr;
1827 sum = c[0] + c[1] * xn;
1828 eps = gnm_abs (GNM_EPSILON / 2 * sum);
1829 for (ui = 2; ui < G_N_ELEMENTS (c); ui++) {
1830 gnm_float t;
1831 xn *= dx;
1832 t = c[ui] * xn;
1833 sum += t;
1834 if (gnm_abs (t) < eps)
1835 break;
1836 }
1837
1838 return sum / x;
1839 }
1840
1841 static gnm_float
gnm_digamma_asymp(gnm_float x)1842 gnm_digamma_asymp (gnm_float x)
1843 {
1844 // Use asympototic series for exp(digamma(x+1/2))
1845 // The use of exp here makes for less cancellation. Note, that the
1846 // asympototic series for plain digamma has a log(x) term, so we
1847 // need the log call anyway. The use of +1/2 makes all the even
1848 // powers go away.
1849
1850 // There are more terms here than will be used in practice
1851 static const gnm_float c[] = {
1852 1, // x
1853 GNM_const(+0.04166666666666666666666667), // x^-1
1854 GNM_const(-0.006423611111111111111111111), // x^-3
1855 GNM_const(+0.003552482914462081128747795), // x^-5
1856 GNM_const(-0.003953557448973030570252792), // x^-7
1857 GNM_const(+0.007365033269308668975914346), // x^-9
1858 GNM_const(-0.02073467582436813806307797), // x^-11
1859 GNM_const(+0.08238185223878776450850024), // x^-13
1860 GNM_const(-0.4396044368600812717750832), // x^-15
1861 GNM_const(+3.034822873180573561262723), // x^-17
1862 GNM_const(-26.32566091447594628148156) // x^-19
1863 };
1864
1865 gnm_float z = x - 0.5, zm2 = 1 / (z * z), zn = z;
1866 gnm_float eps = GNM_EPSILON * z;
1867 gnm_float sum = z;
1868 unsigned ui;
1869
1870 for (ui = 1; ui < G_N_ELEMENTS (c); ui++) {
1871 gnm_float t;
1872 zn *= zm2;
1873 t = c[ui] * zn;
1874 sum += t;
1875 if (gnm_abs (t) < eps)
1876 break;
1877 }
1878
1879 return gnm_log (sum);
1880 }
1881
1882
1883 /**
1884 * gnm_digamma:
1885 * @x: a number
1886 *
1887 * Returns: the digamma function evaluated at @x.
1888 */
1889 gnm_float
gnm_digamma(gnm_float x)1890 gnm_digamma (gnm_float x)
1891 {
1892 // x0 = x0a + x0b is the positive root
1893 gnm_float x0a = GNM_const(1.4616321449683622457627052426687441766262054443359375);
1894 gnm_float x0b = GNM_const(9.549995429965697715184199075967050885129598840859878644035380181024307499273372559036557380022743e-17);
1895
1896 if (gnm_isnan (x))
1897 return x;
1898
1899 if (x <= 0) {
1900 if (x == gnm_floor (x))
1901 return gnm_nan; // Including infinite
1902
1903 // Reflection. Not ideal near zeros
1904 return gnm_digamma (1 - x) - M_PIgnum * gnm_cotpi (x);
1905 }
1906
1907 if (x < x0a - 1)
1908 // Single step up. No cancellation as digamma is negative
1909 // at x+1.
1910 return gnm_digamma (x + 1) - 1 / x;
1911
1912 if (x < x0a - 1.0 / 3.0)
1913 // Series for range [0.46;1.13]
1914 return gnm_digamma_series_1 (x);
1915
1916 if (x < x0a + 1.0 / 3.0)
1917 // Series for range [1.13;1.79] around x0
1918 // Take extra care to compute the difference to x0 with a high-
1919 // precision version of x0
1920 return gnm_digamma_series_2 (x, (x - x0a) - x0b);
1921
1922 if (x < x0a + 1)
1923 // Series for range [1.79;2.46]
1924 return gnm_digamma_series_3 (x);
1925
1926 if (x < 20) {
1927 // Step down to series 2 or 3. All terms are positive so no
1928 // cancellation.
1929 gnm_float sum = 0;
1930 while (x > x0a + 1) {
1931 x--;
1932 sum += 1.0 / x;
1933 }
1934 return sum + gnm_digamma (x);
1935 }
1936
1937 return gnm_digamma_asymp (x);
1938 }
1939
1940 /* ------------------------------------------------------------------------- */
1941