single-gamic.c
1 //! @file single-gamic.c
2 //! @author J. Marcel van der Veer
3 //!
4 //! @section Copyright
5 //!
6 //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
7 //! Copyright 2001-2023 J. Marcel van der Veer [algol68g@xs4all.nl].
8 //!
9 //! @section License
10 //!
11 //! This program is free software; you can redistribute it and/or modify it
12 //! under the terms of the GNU General Public License as published by the
13 //! Free Software Foundation; either version 3 of the License, or
14 //! (at your option) any later version.
15 //!
16 //! This program is distributed in the hope that it will be useful, but
17 //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
19 //! more details. You should have received a copy of the GNU General Public
20 //! License along with this program. If not, see [http://www.gnu.org/licenses/].
21
22 //! @section Synopsis
23 //!
24 //! REAL generalised incomplete gamma function.
25
26 // Generalised incomplete gamma code in this file was downloaded from
27 // [http://helios.mi.parisdescartes.fr/~rabergel/]
28 // and adapted for Algol 68 Genie.
29 //
30 // Reference:
31 // Rémy Abergel, Lionel Moisan. Fast and accurate evaluation of a
32 // generalized incomplete gamma function. 2019. hal-01329669v2
33 //
34 // Original source code copyright and license:
35 //
36 // DELTAGAMMAINC Fast and Accurate Evaluation of a Generalized Incomplete Gamma
37 // Function. Copyright (C) 2016 Remy Abergel (remy.abergel AT gmail.com), Lionel
38 // Moisan (Lionel.Moisan AT parisdescartes.fr).
39 //
40 // This file is a part of the DELTAGAMMAINC software, dedicated to the
41 // computation of a generalized incomplete gammafunction. See the Companion paper
42 // for a complete description of the algorithm.
43 //
44 // ``Fast and accurate evaluation of a generalized incomplete gamma function''
45 // (Rémy Abergel, Lionel Moisan), preprint MAP5 nº2016-14, revision 1.
46 //
47 // This program is free software: you can redistribute it and/or modify it under
48 // the terms of the GNU General Public License as published by the Free Software
49 // Foundation, either version 3 of the License, or (at your option) any later
50 // version.
51 //
52 // This program is distributed in the hope that it will be useful, but WITHOUT
53 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
54 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
55 // details.
56 //
57 // You should have received a copy of the GNU General Public License along with
58 // this program. If not, see [http://www.gnu.org/licenses/].
59
60 // References
61 //
62 // R. Abergel and L. Moisan. 2016. Fast and accurate evaluation of a
63 // generalized incomplete gamma function, preprint MAP5 nº2016-14, revision 1
64 //
65 // Rémy Abergel, Lionel Moisan. Fast and accurate evaluation of a
66 // generalized incomplete gamma function. 2019. hal-01329669v2
67 //
68 // F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark
69 // (Eds.). 2010. NIST Handbook of Mathematical Functions. Cambridge University
70 // Press. (see online version at [http://dlmf.nist.gov/])
71 //
72 // W. H. Press, S. A. Teukolsky, W. T. Vetterling, and
73 // B. P. Flannery. 1992. Numerical recipes in C: the art of scientific
74 // computing (2nd ed.).
75 //
76 // G. R. Pugh, 2004. An analysis of the Lanczos Gamma approximation (phd
77 // thesis)
78
79 #include "a68g.h"
80 #include "a68g-genie.h"
81 #include "a68g-prelude.h"
82 #include "a68g-lib.h"
83 #include "a68g-double.h"
84 #include "a68g-mp.h"
85
86 #define ITMAX 1000000000 // Maximum allowed number of iterations
87 #define DPMIN DBL_MIN // Number near the smallest representable double-point number
88 #define EPS DBL_EPSILON // Machine epsilon
89 #define NITERMAX_ROMBERG 15 // Maximum allowed number of Romberg iterations
90 #define TOL_ROMBERG 0.1 // Tolerance factor used to stop the Romberg iterations
91 #define TOL_DIFF 0.2 // Tolerance factor used for the approximation of I_{x,y}^{mu,p} using differences
92
93 // plim: compute plim (x), the limit of the partition of the domain (p,x)
94 // detailed in the paper.
95 //
96 // | x if 0 < x
97 // |
98 // plim (x) = < 0 if -9 <= x <= 0
99 // |
100 // | 5.*sqrt (|x|)-5. otherwise
101
102 REAL_T plim (REAL_T x)
103 {
104 return (x >= 0) ? x : ((x >= -9) ? 0 : 5 * sqrt (-x) - 5);
105 }
106
107 //! @brief compute G(p,x) in the domain x <= p using a continued fraction
108 //
109 // p >= 0
110 // x <= p
111
112 void G_cfrac_lower (REAL_T * Gcfrac, REAL_T p, REAL_T x)
113 {
114 REAL_T c, d, del, f, an, bn;
115 INT_T k, n;
116 // deal with special case
117 if (x == 0) {
118 *Gcfrac = 0;
119 return;
120 }
121 // Evaluate the continued fraction using Modified Lentz's method. However,
122 // as detailed in the paper, perform manually the first pass (n=1), of the
123 // initial Modified Lentz's method.
124 an = 1;
125 bn = p;
126 f = an / bn;
127 c = an / DPMIN;
128 d = 1 / bn;
129 n = 2;
130 do {
131 k = n / 2;
132 an = (n & 1 ? k : -(p - 1 + k)) * x;
133 bn++;
134 d = an * d + bn;
135 if (d == 0) {
136 d = DPMIN;
137 }
138 c = bn + an / c;
139 if (c == 0) {
140 c = DPMIN;
141 }
142 d = 1 / d;
143 del = d * c;
144 f *= del;
145 n++;
146 }
147 while ((a68_abs (del - 1.0) >= EPS) && (n < ITMAX));
148 *Gcfrac = f;
149 }
150
151 //! @brief compute the G-function in the domain x < 0 and |x| < max (1,p-1)
152 // using a recursive integration by parts relation.
153 // This function cannot be used when mu > 0.
154 //
155 // p > 0, integer
156 // x < 0, |x| < max (1,p-1)
157
158 void G_ibp (REAL_T * Gibp, REAL_T p, REAL_T x)
159 {
160 REAL_T t, tt, c, d, s, del;
161 INT_T l;
162 BOOL_T odd, stop;
163 t = a68_abs (x);
164 tt = 1 / (t * t);
165 odd = (INT_T) a68_int (p) % 2 != 0;
166 c = 1 / t;
167 d = (p - 1);
168 s = c * (t - d);
169 l = 0;
170 do {
171 c *= d * (d - 1) * tt;
172 d -= 2;
173 del = c * (t - d);
174 s += del;
175 l++;
176 stop = a68_abs (del) < a68_abs (s) * EPS;
177 }
178 while ((l < floor ((p - 2) / 2)) && !stop);
179 if (odd && !stop) {
180 s += d * c / t;
181 }
182 *Gibp = ((odd ? -1 : 1) * a68_exp (-t + lgamma (p) - (p - 1) * a68_ln (t)) + s) / t;
183 }
184
185 //! @brief compute the G-function in the domain x > p using a
186 // continued fraction.
187 //
188 // p > 0
189 // x > p, or x = +infinity
190
191 void G_cfrac_upper (REAL_T * Gcfrac, REAL_T p, REAL_T x)
192 {
193 // Special case
194 if (a68_isinf (x)) {
195 *Gcfrac = 0;
196 return;
197 }
198 // Evaluate the continued fraction using Modified Lentz's method. However,
199 // as detailed in the paper, perform manually the first pass (n=1), of the
200 // initial Modified Lentz's method.
201 REAL_T an = 1, bn = x + 1 - p;
202 REAL_T c, d, del, f;
203 BOOL_T t = (bn != 0);
204 INT_T i, n;
205 if (t) {
206 // b{1} is non-zero
207 f = an / bn;
208 c = an / DPMIN;
209 d = 1 / bn;
210 n = 2;
211 } else {
212 // b{1}=0 but b{2} is non-zero, compute Mcfrac = a{1}/f with f = a{2}/(b{2}+) a{3}/(b{3}+) ...
213 an = -(1 - p);
214 bn = x + 3 - p;
215 f = an / bn;
216 c = an / DPMIN;
217 d = 1 / bn;
218 n = 3;
219 }
220 i = n - 1;
221 do {
222 an = -i * (i - p);
223 bn += 2;
224 d = an * d + bn;
225 if (d == 0) {
226 d = DPMIN;
227 }
228 c = bn + an / c;
229 if (c == 0) {
230 c = DPMIN;
231 }
232 d = 1 / d;
233 del = d * c;
234 f *= del;
235 i++;
236 n++;
237 }
238 while ((a68_abs (del - 1.0) >= EPS) && (n < ITMAX));
239 *Gcfrac = t ? f : 1 / f;
240 }
241
242 //! @brief compute G : (p,x) --> R defined as follows
243 //
244 // if x <= p:
245 // G(p,x) = exp (x-p*ln (|x|)) * integral of s^{p-1} * exp (-sign (x)*s) ds from s = 0 to |x|
246 // otherwise:
247 // G(p,x) = exp (x-p*ln (|x|)) * integral of s^{p-1} * exp (-s) ds from s = x to infinity
248 //
249 // p > 0
250 // x is a real number or +infinity.
251
252 void G_func (REAL_T * G, REAL_T p, REAL_T x)
253 {
254 if (p >= plim (x)) {
255 G_cfrac_lower (G, p, x);
256 } else if (x < 0) {
257 G_ibp (G, p, x);
258 } else {
259 G_cfrac_upper (G, p, x);
260 }
261 }
262
263 //! @brief iteration of the Romberg approximation of I_{x,y}^{mu,p}
264
265 void romberg_iterations (REAL_T * R, REAL_T sigma, INT_T n, REAL_T x, REAL_T y, REAL_T mu, REAL_T p, REAL_T h, REAL_T pow2)
266 {
267 INT_T j, m;
268 REAL_T sum, xx;
269 INT_T adr0_prev = ((n - 1) * n) / 2;
270 INT_T adr0 = (n * (n + 1)) / 2;
271 for (sum = 0, j = 1; j <= pow2; j++) {
272 xx = x + ((y - x) * (2 * j - 1)) / (2 * pow2);
273 sum += a68_exp (-mu * xx + (p - 1) * a68_ln (xx) - sigma);
274 }
275 R[adr0] = 0.5 * R[adr0_prev] + h * sum;
276 REAL_T pow4 = 4;
277 for (m = 1; m <= n; m++) {
278 R[adr0 + m] = (pow4 * R[adr0 + (m - 1)] - R[adr0_prev + (m - 1)]) / (pow4 - 1);
279 pow4 *= 4;
280 }
281 }
282
283 //! @ compute I_{x,y}^{mu,p} using a Romberg approximation.
284 // Compute rho and sigma so I_{x,y}^{mu,p} = rho * exp (sigma)
285
286 void romberg_estimate (REAL_T * rho, REAL_T * sigma, REAL_T x, REAL_T y, REAL_T mu, REAL_T p)
287 {
288 REAL_T *R = (REAL_T *) get_heap_space (((NITERMAX_ROMBERG + 1) * (NITERMAX_ROMBERG + 2)) / 2 * sizeof (REAL_T));
289 ASSERT (R != NULL);
290 // Initialization (n=1)
291 *sigma = -mu * y + (p - 1) * a68_ln (y);
292 R[0] = 0.5 * (y - x) * (a68_exp (-mu * x + (p - 1) * a68_ln (x) - (*sigma)) + 1);
293 // Loop for n > 0
294 REAL_T relneeded = EPS / TOL_ROMBERG;
295 INT_T adr0 = 0;
296 INT_T n = 1;
297 REAL_T h = (y - x) / 2; // n=1, h = (y-x)/2^n
298 REAL_T pow2 = 1; // n=1; pow2 = 2^(n-1)
299 if (NITERMAX_ROMBERG >= 1) {
300 REAL_T relerr;
301 do {
302 romberg_iterations (R, *sigma, n, x, y, mu, p, h, pow2);
303 h /= 2;
304 pow2 *= 2;
305 adr0 = (n * (n + 1)) / 2;
306 relerr = a68_abs ((R[adr0 + n] - R[adr0 + n - 1]) / R[adr0 + n]);
307 n++;
308 } while (n <= NITERMAX_ROMBERG && relerr > relneeded);
309 }
310 // save Romberg estimate and free memory
311 *rho = R[adr0 + (n - 1)];
312 a68_free (R);
313 }
314
315 //! @brief compute generalized incomplete gamma function I_{x,y}^{mu,p}
316 //
317 // I_{x,y}^{mu,p} = integral from x to y of s^{p-1} * exp (-mu*s) ds
318 //
319 // This procedure computes (rho, sigma) described below.
320 // The approximated value of I_{x,y}^{mu,p} is I = rho * exp (sigma)
321 //
322 // mu is a real number non equal to zero
323 // (in general we take mu = 1 or -1 but any nonzero real number is allowed)
324 //
325 // x, y are two numbers with 0 <= x <= y <= +infinity,
326 // (the setting y=+infinity is allowed only when mu > 0)
327 //
328 // p is a real number > 0, p must be an integer when mu < 0.
329
330 void deltagammainc (REAL_T * rho, REAL_T * sigma, REAL_T x, REAL_T y, REAL_T mu, REAL_T p)
331 {
332 // Particular cases
333 if (a68_isinf (x) && a68_isinf (y)) {
334 *rho = 0;
335 *sigma = a68_neginf ();
336 return;
337 } else if (x == y) {
338 *rho = 0;
339 *sigma = a68_neginf ();
340 return;
341 }
342 if (x == 0 && a68_isinf (y)) {
343 *rho = 1;
344 (*sigma) = lgamma (p) - p * a68_ln (mu);
345 return;
346 }
347 // Initialization
348 REAL_T mA, mB, mx, my, nA, nB, nx, ny;
349 G_func (&mx, p, mu * x);
350 nx = (a68_isinf (x) ? a68_neginf () : -mu * x + p * a68_ln (x));
351 G_func (&my, p, mu * y);
352 ny = (a68_isinf (y) ? a68_neginf () : -mu * y + p * a68_ln (y));
353 // Compute (mA,nA) and (mB,nB) such as I_{x,y}^{mu,p} can be
354 // approximated by the difference A-B, where A >= B >= 0, A = mA*exp (nA) an
355 // B = mB*exp (nB). When the difference involves more than one digit loss due to
356 // cancellation errors, the integral I_{x,y}^{mu,p} is evaluated using the
357 // Romberg approximation method.
358
359 if (mu < 0) {
360 mA = my;
361 nA = ny;
362 mB = mx;
363 nB = nx;
364 } else {
365 if (p < plim (mu * x)) {
366 mA = mx;
367 nA = nx;
368 mB = my;
369 nB = ny;
370 } else if (p < plim (mu * y)) {
371 mA = 1;
372 nA = lgamma (p) - p * a68_ln (mu);
373 nB = fmax (nx, ny);
374 mB = mx * a68_exp (nx - nB) + my * a68_exp (ny - nB);
375 } else {
376 mA = my;
377 nA = ny;
378 mB = mx;
379 nB = nx;
380 }
381 }
382 // Compute (rho,sigma) such that rho*exp (sigma) = A-B
383 *rho = mA - mB * a68_exp (nB - nA);
384 *sigma = nA;
385 // If the difference involved a significant loss of precision, compute Romberg estimate.
386 if (!a68_isinf (y) && ((*rho) / mA < TOL_DIFF)) {
387 romberg_estimate (rho, sigma, x, y, mu, p);
388 }
389 }
390
391 // A68G Driver routines
392
393 //! @brief PROC gamma inc g = (REAL p, x, y, mu) REAL
394
395 void genie_gamma_inc_g_real (NODE_T * n)
396 {
397 A68_REAL x, y, mu, p;
398 POP_OBJECT (n, &mu, A68_REAL);
399 POP_OBJECT (n, &y, A68_REAL);
400 POP_OBJECT (n, &x, A68_REAL);
401 POP_OBJECT (n, &p, A68_REAL);
402 REAL_T rho, sigma;
403 deltagammainc (&rho, &sigma, VALUE (&x), VALUE (&y), VALUE (&mu), VALUE (&p));
404 PUSH_VALUE (n, rho * a68_exp (sigma), A68_REAL);
405 }
406
407 //! @brief PROC gamma inc f = (REAL p, x) REAL
408
409 void genie_gamma_inc_f_real (NODE_T * n)
410 {
411 A68_REAL x, p;
412 POP_OBJECT (n, &x, A68_REAL);
413 POP_OBJECT (n, &p, A68_REAL);
414 REAL_T rho, sigma;
415 deltagammainc (&rho, &sigma, VALUE (&x), a68_posinf (), 1, VALUE (&p));
416 PUSH_VALUE (n, rho * a68_exp (sigma), A68_REAL);
417 }
418
419 //! @brief PROC gamma inc = (REAL p, x) REAL
420
421 void genie_gamma_inc_h_real (NODE_T * n)
422 {
423 #if (A68_LEVEL >= 3) && defined (HAVE_GNU_MPFR)
424 genie_gamma_inc_real_mpfr (n);
425 #else
426 genie_gamma_inc_f_real (n);
427 #endif
428 }
429
430 //! @brief PROC gamma inc gf = (REAL p, x) REAL
431
432 void genie_gamma_inc_gf_real (NODE_T * q)
433 {
434 // if x <= p: G(p,x) = exp (x-p*ln (|x|)) * integral over [0,|x|] of s^{p-1} * exp (-sign (x)*s) ds
435 // otherwise: G(p,x) = exp (x-p*ln (x)) * integral over [x,inf] of s^{p-1} * exp (-s) ds
436 A68_REAL x, p;
437 POP_OBJECT (q, &x, A68_REAL);
438 POP_OBJECT (q, &p, A68_REAL);
439 REAL_T G;
440 G_func (&G, VALUE (&p), VALUE (&x));
441 PUSH_VALUE (q, G, A68_REAL);
442 }
