double-gamic.c
1 //! @file double-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 //! LONG 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
81 #if (A68_LEVEL >= 3)
82
83 #include "a68g-double.h"
84 #include "a68g-genie.h"
85 #include "a68g-lib.h"
86 #include "a68g-mp.h"
87 #include "a68g-prelude.h"
88 #include "a68g-quad.h"
89
90 #define DPMIN FLT128_MIN // Number near the smallest representable double-point number
91 #define EPS FLT128_EPSILON // Machine epsilon
92 #define ITMAX 1000000000 // Maximum allowed number of iterations
93 #define NITERMAX_ROMBERG 15 // Maximum allowed number of Romberg iterations
94 #define TOL_DIFF 0.2q // Tolerance factor used for the approximation of I_{x,y}^{mu,p} using differences
95 #define TOL_ROMBERG 0.1q // Tolerance factor used to stop the Romberg iterations
96
97 // double_plim: compute plim (x), the limit of the partition of the domain (p,x)
98 // detailed in the paper.
99 //
100 // | x if 0 < x
101 // |
102 // plim (x) = < 0 if -9 <= x <= 0
103 // |
104 // | 5.*sqrt (|x|)-5. otherwise
105
106 DOUBLE_T double_plim (DOUBLE_T x)
107 {
108 return (x >= 0.0q) ? x : ((x >= -9.0q) ? 0.0q : 5.0q * sqrt (-x) - 5.0q);
109 }
110
111 //! @brief compute G(p,x) in the domain x <= p using a continued fraction
112 //
113 // p >= 0
114 // x <= p
115
116 void double_G_cfrac_lower (DOUBLE_T * Gcfrac, DOUBLE_T p, DOUBLE_T x)
117 {
118 // deal with special case
119 if (x == 0.0q) {
120 *Gcfrac = 0.0q;
121 return;
122 }
123 // Evaluate the continued fraction using Modified Lentz's method. However,
124 // as detailed in the paper, perform manually the first pass (n=1), of the
125 // initial Modified Lentz's method.
126 DOUBLE_T an = 1.0q, bn = p, del;
127 DOUBLE_T f = an / bn, c = an / DPMIN, d = 1.0q / bn;
128 INT_T n = 2;
129 do {
130 INT_T k = n / 2;
131 an = (n & 1 ? k : -(p - 1.0q + k)) * x;
132 bn++;
133 d = an * d + bn;
134 if (d == 0.0q) {
135 d = DPMIN;
136 }
137 c = bn + an / c;
138 if (c == 0.0q) {
139 c = DPMIN;
140 }
141 d = 1.0q / d;
142 del = d * c;
143 f *= del;
144 n++;
145 }
146 while ((fabsq (del - 1.0q) >= EPS) && (n < ITMAX));
147 *Gcfrac = f;
148 }
149
150 //! @brief compute the G-function in the domain x < 0 and |x| < max (1,p-1)
151 // using a recursive integration by parts relation.
152 // This function cannot be used when mu > 0.
153 //
154 // p > 0, integer
155 // x < 0, |x| < max (1,p-1)
156
157 void double_G_ibp (DOUBLE_T * Gibp, DOUBLE_T p, DOUBLE_T x)
158 {
159 BOOL_T odd = (INT_T) (p) % 2 != 0;
160 DOUBLE_T t = fabsq (x), del;
161 DOUBLE_T tt = 1.0q / (t * t), c = 1.0q / t, d = (p - 1.0q);
162 DOUBLE_T s = c * (t - d);
163 INT_T l = 0;
164 BOOL_T stop;
165 do {
166 c *= d * (d - 1.0q) * tt;
167 d -= 2.0q;
168 del = c * (t - d);
169 s += del;
170 l++;
171 stop = fabsq (del) < fabsq (s) * EPS;
172 }
173 while ((l < floorq ((p - 2.0q) / 2.0q)) && !stop);
174 if (odd && !stop) {
175 s += d * c / t;
176 }
177 *Gibp = ((odd ? -1.0q : 1.0q) * expq (-t + lgammaq (p) - (p - 1.0q) * logq (t)) + s) / t;
178 }
179
180 //! @brief compute the G-function in the domain x > p using a
181 // continued fraction.
182 //
183 // p > 0
184 // x > p, or x = +infinity
185
186 void double_G_cfrac_upper (DOUBLE_T * Gcfrac, DOUBLE_T p, DOUBLE_T x)
187 {
188 DOUBLE_T c, d, del, f, an, bn;
189 INT_T i, n;
190 // Special case
191 if (isinfq (x)) {
192 *Gcfrac = 0.0q;
193 return;
194 }
195 // Evaluate the continued fraction using Modified Lentz's method. However,
196 // as detailed in the paper, perform manually the first pass (n=1), of the
197 // initial Modified Lentz's method.
198 an = 1.0q;
199 bn = x + 1.0q - p;
200 BOOL_T t = (bn != 0.0q);
201 if (t) {
202 // b{1} is non-zero
203 f = an / bn;
204 c = an / DPMIN;
205 d = 1.0q / bn;
206 n = 2;
207 } else {
208 // b{1}=0 but b{2} is non-zero, compute Mcfrac = a{1}/f with f = a{2}/(b{2}+) a{3}/(b{3}+) ...
209 an = -(1.0q - p);
210 bn = x + 3.0q - p;
211 f = an / bn;
212 c = an / DPMIN;
213 d = 1.0q / bn;
214 n = 3;
215 }
216 i = n - 1;
217 do {
218 an = -i * (i - p);
219 bn += 2.0q;
220 d = an * d + bn;
221 if (d == 0.0q) {
222 d = DPMIN;
223 }
224 c = bn + an / c;
225 if (c == 0.0q) {
226 c = DPMIN;
227 }
228 d = 1.0q / d;
229 del = d * c;
230 f *= del;
231 i++;
232 n++;
233 }
234 while ((fabsq (del - 1.0q) >= EPS) && (n < ITMAX));
235 *Gcfrac = t ? f : 1.0q / f;
236 }
237
238 //! @brief compute G : (p,x) --> R defined as follows
239 //
240 // if x <= p:
241 // G(p,x) = exp (x-p*ln (|x|)) * integral of s^{p-1} * exp (-sign (x)*s) ds from s = 0 to |x|
242 // otherwise:
243 // G(p,x) = exp (x-p*ln (|x|)) * integral of s^{p-1} * exp (-s) ds from s = x to infinity
244 //
245 // p > 0
246 // x is a real number or +infinity.
247
248 void double_G_func (DOUBLE_T * G, DOUBLE_T p, DOUBLE_T x)
249 {
250 if (p >= double_plim (x)) {
251 double_G_cfrac_lower (G, p, x);
252 } else if (x < 0.0q) {
253 double_G_ibp (G, p, x);
254 } else {
255 double_G_cfrac_upper (G, p, x);
256 }
257 }
258
259 //! @brief iteration of the Romberg approximation of I_{x,y}^{mu,p}
260
261 void double_romberg_iterations (DOUBLE_T * R, DOUBLE_T sigma, INT_T n, DOUBLE_T x, DOUBLE_T y, DOUBLE_T mu, DOUBLE_T p, DOUBLE_T h, DOUBLE_T pow2)
262 {
263 INT_T adr0_prev = ((n - 1) * n) / 2, adr0 = (n * (n + 1)) / 2;
264 QUAD_T sum = QUAD_REAL_ZERO;
265 for (INT_T j = 1; j <= pow2; j++) {
266 DOUBLE_T xx = x + ((y - x) * (2.0q * j - 1.0q)) / (2.0q * pow2);
267 QUAD_T f = double_real_to_quad_real (expq (-mu * xx + (p - 1.0q) * logq (xx) - sigma));
268 sum = _add_quad_real_ (sum, f);
269 }
270 R[adr0] = 0.5q * R[adr0_prev] + h * quad_real_to_double_real (sum);
271 DOUBLE_T pow4 = 4.0q;
272 for (INT_T m = 1; m <= n; m++) {
273 R[adr0 + m] = (pow4 * R[adr0 + (m - 1)] - R[adr0_prev + (m - 1)]) / (pow4 - 1.0q);
274 pow4 *= 4.0q;
275 }
276 }
277
278 //! @ compute I_{x,y}^{mu,p} using a Romberg approximation.
279 // Compute rho and sigma so I_{x,y}^{mu,p} = rho * exp (sigma)
280
281 void double_romberg_estimate (DOUBLE_T * rho, DOUBLE_T * sigma, DOUBLE_T x, DOUBLE_T y, DOUBLE_T mu, DOUBLE_T p)
282 {
283 DOUBLE_T *R = (DOUBLE_T *) get_heap_space (((NITERMAX_ROMBERG + 1) * (NITERMAX_ROMBERG + 2)) / 2 * sizeof (DOUBLE_T));
284 ASSERT (R != NULL);
285 // Initialization (n=1)
286 *sigma = -mu * y + (p - 1.0q) * logq (y);
287 R[0] = 0.5q * (y - x) * (expq (-mu * x + (p - 1.0q) * logq (x) - (*sigma)) + 1.0q);
288 // Loop for n > 0
289 DOUBLE_T relneeded = EPS / TOL_ROMBERG;
290 INT_T adr0 = 0;
291 INT_T n = 1;
292 DOUBLE_T h = (y - x) / 2.0q; // n=1, h = (y-x)/2^n
293 DOUBLE_T pow2 = 1.0q; // n=1; pow2 = 2^(n-1)
294 if (NITERMAX_ROMBERG >= 1) {
295 DOUBLE_T relerr;
296 do {
297 double_romberg_iterations (R, *sigma, n, x, y, mu, p, h, pow2);
298 h /= 2.0q;
299 pow2 *= 2.0q;
300 adr0 = (n * (n + 1)) / 2;
301 relerr = fabsq ((R[adr0 + n] - R[adr0 + n - 1]) / R[adr0 + n]);
302 n++;
303 } while (n <= NITERMAX_ROMBERG && relerr > relneeded);
304 }
305 // save Romberg estimate and free memory
306 *rho = R[adr0 + (n - 1)];
307 a68_free (R);
308 }
309
310 //! @brief compute generalized incomplete gamma function I_{x,y}^{mu,p}
311 //
312 // I_{x,y}^{mu,p} = integral from x to y of s^{p-1} * exp (-mu*s) ds
313 //
314 // This procedure computes (rho, sigma) described below.
315 // The approximated value of I_{x,y}^{mu,p} is I = rho * exp (sigma)
316 //
317 // mu is a real number non equal to zero
318 // (in general we take mu = 1 or -1 but any nonzero real number is allowed)
319 //
320 // x, y are two numbers with 0 <= x <= y <= +infinity,
321 // (the setting y=+infinity is allowed only when mu > 0)
322 //
323 // p is a real number > 0, p must be an integer when mu < 0.
324
325 void deltagammainc_double_real (DOUBLE_T * rho, DOUBLE_T * sigma, DOUBLE_T x, DOUBLE_T y, DOUBLE_T mu, DOUBLE_T p)
326 {
327 DOUBLE_T mA, mB, mx, my, nA, nB, nx, ny;
328 // Particular cases
329 if (isinfq (x) && isinfq (y)) {
330 *rho = 0.0q;
331 *sigma = a68_dneginf ();
332 return;
333 } else if (x == y) {
334 *rho = 0.0q;
335 *sigma = a68_dneginf ();
336 return;
337 }
338 if (x == 0.0q && isinfq (y)) {
339 *rho = 1.0q;
340 (*sigma) = lgammaq (p) - p * logq (mu);
341 return;
342 }
343 // Initialization
344 double_G_func (&mx, p, mu * x);
345 nx = (isinfq (x) ? a68_dneginf () : -mu * x + p * logq (x));
346 double_G_func (&my, p, mu * y);
347 ny = (isinfq (y) ? a68_dneginf () : -mu * y + p * logq (y));
348
349 // Compute (mA,nA) and (mB,nB) such as I_{x,y}^{mu,p} can be
350 // approximated by the difference A-B, where A >= B >= 0, A = mA*exp (nA) an
351 // B = mB*exp (nB). When the difference involves more than one digit loss due to
352 // cancellation errors, the integral I_{x,y}^{mu,p} is evaluated using the
353 // Romberg approximation method.
354
355 if (mu < 0.0q) {
356 mA = my;
357 nA = ny;
358 mB = mx;
359 nB = nx;
360 } else {
361 if (p < double_plim (mu * x)) {
362 mA = mx;
363 nA = nx;
364 mB = my;
365 nB = ny;
366 } else if (p < double_plim (mu * y)) {
367 mA = 1.0q;
368 nA = lgammaq (p) - p * logq (mu);
369 nB = fmax (nx, ny);
370 mB = mx * expq (nx - nB) + my * expq (ny - nB);
371 } else {
372 mA = my;
373 nA = ny;
374 mB = mx;
375 nB = nx;
376 }
377 }
378 // Compute (rho,sigma) such that rho*exp (sigma) = A-B
379 *rho = mA - mB * expq (nB - nA);
380 *sigma = nA;
381 // If the difference involved a significant loss of precision, compute Romberg estimate.
382 if (!isinfq (y) && ((*rho) / mA < TOL_DIFF)) {
383 double_romberg_estimate (rho, sigma, x, y, mu, p);
384 }
385 }
386
387 // A68G Driver routines
388
389 //! @brief PROC long gamma inc g = (LONG REAL p, x, y, mu) LONG REAL
390
391 void genie_gamma_inc_g_double_real (NODE_T * n)
392 {
393 A68_LONG_REAL x, y, mu, p;
394 POP_OBJECT (n, &mu, A68_LONG_REAL);
395 POP_OBJECT (n, &y, A68_LONG_REAL);
396 POP_OBJECT (n, &x, A68_LONG_REAL);
397 POP_OBJECT (n, &p, A68_LONG_REAL);
398 DOUBLE_T rho, sigma;
399 deltagammainc_double_real (&rho, &sigma, VALUE (&x).f, VALUE (&y).f, VALUE (&mu).f, VALUE (&p).f);
400 PUSH_VALUE (n, dble (rho * expq (sigma)), A68_LONG_REAL);
401 }
402
403 //! @brief PROC long gamma inc f = (LONG REAL p, x) LONG REAL
404
405 void genie_gamma_inc_f_double_real (NODE_T * n)
406 {
407 A68_LONG_REAL x, p;
408 POP_OBJECT (n, &x, A68_LONG_REAL);
409 POP_OBJECT (n, &p, A68_LONG_REAL);
410 DOUBLE_T rho, sigma;
411 deltagammainc_double_real (&rho, &sigma, VALUE (&x).f, a68_dposinf (), 1.0q, VALUE (&p).f);
412 PUSH_VALUE (n, dble (rho * expq (sigma)), A68_LONG_REAL);
413 }
414
415 //! @brief PROC long gamma inc gf = (LONG REAL p, x) LONG REAL
416
417 void genie_gamma_inc_gf_double_real (NODE_T * q)
418 {
419 // if x <= p: G(p,x) = exp (x-p*ln (|x|)) * integral over [0,|x|] of s^{p-1} * exp (-sign (x)*s) ds
420 // otherwise: G(p,x) = exp (x-p*ln (x)) * integral over [x,inf] of s^{p-1} * exp (-s) ds
421 A68_LONG_REAL x, p;
422 POP_OBJECT (q, &x, A68_LONG_REAL);
423 POP_OBJECT (q, &p, A68_LONG_REAL);
424 DOUBLE_T G;
425 double_G_func (&G, VALUE (&p).f, VALUE (&x).f);
426 PUSH_VALUE (q, dble (G), A68_LONG_REAL);
427 }
428
429 //! @brief PROC long gamma inc = (LONG REAL p, x) LONG REAL
430
431 void genie_gamma_inc_h_double_real (NODE_T * n)
432 {
433 #if (A68_LEVEL >= 3) && defined (HAVE_GNU_MPFR)
434 genie_gamma_inc_double_real_mpfr (n);
435 #else
436 genie_gamma_inc_f_double_real (n);
437 #endif
438 }
439
440 #endif
