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