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-2024 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 //
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 // DELTAGAMMAINC Fast and Accurate Evaluation of a Generalized Incomplete Gamma
38 // Function. Copyright (C) 2016 Remy Abergel (remy.abergel AT gmail.com), Lionel
39 // Moisan (Lionel.Moisan AT parisdescartes.fr).
40 //
41 // This file is a part of the DELTAGAMMAINC software, dedicated to the
42 // computation of a generalized incomplete gammafunction. See the Companion paper
43 // for a complete description of the algorithm.
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 // This program is distributed in the hope that it will be useful, but WITHOUT
52 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
53 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
54 // details.
55 //
56 // You should have received a copy of the GNU General Public License along with
57 // this program. If not, see [http://www.gnu.org/licenses/].
58
59 // References
60 // R. Abergel and L. Moisan. 2016. Fast and accurate evaluation of a
61 // generalized incomplete gamma function, preprint MAP5 nº2016-14, revision 1
62 // Rémy Abergel, Lionel Moisan. Fast and accurate evaluation of a
63 // generalized incomplete gamma function. 2019. hal-01329669v2
64 // F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark
65 // (Eds.). 2010. NIST Handbook of Mathematical Functions. Cambridge University
66 // Press. (see online version at [http://dlmf.nist.gov/])
67 // W. H. Press, S. A. Teukolsky, W. T. Vetterling, and
68 // B. P. Flannery. 1992. Numerical recipes in C: the art of scientific
69 // computing (2nd ed.).
70 // G. R. Pugh, 2004. An analysis of the Lanczos Gamma approximation (phd
71 // thesis)
72
73 #include "a68g.h"
74
75 #if (A68_LEVEL >= 3)
76
77 #include "a68g-double.h"
78 #include "a68g-genie.h"
79 #include "a68g-mp.h"
80
81 #define ITMAX 1000000000 // Maximum allowed number of iterations
82 #define NITERMAX_ROMBERG 15 // Maximum allowed number of Romberg iterations
83 #define TOL_DIFF 0.2q // Tolerance factor used for the approximation of I_{x,y}^{mu,p} using differences
84 #define TOL_ROMBERG 0.1q // Tolerance factor used to stop the Romberg iterations
85
86 //! @brief compute G(p,x) in the domain x <= p using a continued fraction
87 // p >= 0
88 // x <= p
89
90 void double_G_cfrac_lower (DOUBLE_T * Gcfrac, DOUBLE_T p, DOUBLE_T x)
91 {
92 // deal with special case
93 if (x == 0.0q) {
94 *Gcfrac = 0.0q;
95 return;
96 }
97 // Evaluate the continued fraction using Modified Lentz's method. However,
98 // as detailed in the paper, perform manually the first pass (n=1), of the
99 // initial Modified Lentz's method.
100 DOUBLE_T an = 1.0q, bn = p, del;
101 DOUBLE_T f = an / bn, c = an / A68_DOUBLE_MIN, d = 1.0q / bn;
102 INT_T n = 2;
103 do {
104 INT_T k = n / 2;
105 an = (n & 1 ? k : -(p - 1.0q + k)) * x;
106 bn++;
107 d = an * d + bn;
108 if (d == 0.0q) {
109 d = A68_DOUBLE_MIN;
110 }
111 c = bn + an / c;
112 if (c == 0.0q) {
113 c = A68_DOUBLE_MIN;
114 }
115 d = 1.0q / d;
116 del = d * c;
117 f *= del;
118 n++;
119 }
120 while ((fabs_double (del - 1.0q) >= A68_DOUBLE_EPS) && (n < ITMAX));
121 *Gcfrac = f;
122 }
123
124 //! @brief compute the G-function in the domain x < 0 and |x| < max (1,p-1)
125 // using a recursive integration by parts relation.
126 // This function cannot be used when mu > 0.
127 // p > 0, integer
128 // x < 0, |x| < max (1,p-1)
129
130 void double_G_ibp (DOUBLE_T * Gibp, DOUBLE_T p, DOUBLE_T x)
131 {
132 BOOL_T odd = (INT_T) (p) % 2 != 0;
133 DOUBLE_T t = fabs_double (x), del;
134 DOUBLE_T tt = 1.0q / (t * t), c = 1.0q / t, d = (p - 1.0q);
135 DOUBLE_T s = c * (t - d);
136 INT_T l = 0;
137 BOOL_T stop;
138 do {
139 c *= d * (d - 1.0q) * tt;
140 d -= 2.0q;
141 del = c * (t - d);
142 s += del;
143 l++;
144 stop = fabs_double (del) < fabs_double (s) * A68_DOUBLE_EPS;
145 }
146 while ((l < floor_double ((p - 2.0q) / 2.0q)) && !stop);
147 if (odd && !stop) {
148 s += d * c / t;
149 }
150 *Gibp = ((odd ? -1.0q : 1.0q) * exp_double (-t + lgamma_double (p) - (p - 1.0q) * log_double (t)) + s) / t;
151 }
152
153 //! @brief compute plim (x), the limit of the partition of the domain (p,x)
154 // detailed in the paper.
155 // | x if 0 < x
156 // |
157 // plim (x) = < 0 if -9 <= x <= 0
158 // |
159 // | 5.*sqrt (|x|)-5. otherwise
160
161 DOUBLE_T double_plim (DOUBLE_T x)
162 {
163 return (x >= 0.0q) ? x : ((x >= -9.0q) ? 0.0q : 5.0q * sqrt (-x) - 5.0q);
164 }
165
166 //! @brief compute the G-function in the domain x > p using a
167 // continued fraction.
168 // p > 0
169 // x > p, or x = +infinity
170
171 void double_G_cfrac_upper (DOUBLE_T * Gcfrac, DOUBLE_T p, DOUBLE_T x)
172 {
173 // Special case
174 if (isinf_double (x)) {
175 *Gcfrac = 0.0q;
176 return;
177 }
178 // Evaluate the continued fraction using Modified Lentz's method. However,
179 // as detailed in the paper, perform manually the first pass (n=1), of the
180 // initial Modified Lentz's method.
181 DOUBLE_T an = 1.0q, bn = x + 1.0q - p, c, d, del, f;
182 BOOL_T t = (bn != 0.0q);
183 INT_T n;
184 if (t) {
185 // b{1} is non-zero
186 f = an / bn;
187 c = an / A68_DOUBLE_MIN;
188 d = 1.0q / bn;
189 n = 2;
190 } else {
191 // b{1}=0 but b{2} is non-zero, compute Mcfrac = a{1}/f with f = a{2}/(b{2}+) a{3}/(b{3}+) ...
192 an = -(1.0q - p);
193 bn = x + 3.0q - p;
194 f = an / bn;
195 c = an / A68_DOUBLE_MIN;
196 d = 1.0q / bn;
197 n = 3;
198 }
199 INT_T i = n - 1;
200 do {
201 an = -i * (i - p);
202 bn += 2.0q;
203 d = an * d + bn;
204 if (d == 0.0q) {
205 d = A68_DOUBLE_MIN;
206 }
207 c = bn + an / c;
208 if (c == 0.0q) {
209 c = A68_DOUBLE_MIN;
210 }
211 d = 1.0q / d;
212 del = d * c;
213 f *= del;
214 i++;
215 n++;
216 }
217 while ((fabs_double (del - 1.0q) >= A68_DOUBLE_EPS) && (n < ITMAX));
218 *Gcfrac = t ? f : 1.0q / f;
219 }
220
221 //! @brief compute G : (p,x) --> R defined as follows
222 // if x <= p:
223 // G(p,x) = exp (x-p*ln (|x|)) * integral of s^{p-1} * exp (-sign (x)*s) ds from s = 0 to |x|
224 // otherwise:
225 // G(p,x) = exp (x-p*ln (|x|)) * integral of s^{p-1} * exp (-s) ds from s = x to infinity
226 // p > 0
227 // x is a real number or +infinity.
228
229 void double_G_func (DOUBLE_T * G, DOUBLE_T p, DOUBLE_T x)
230 {
231 if (p >= double_plim (x)) {
232 double_G_cfrac_lower (G, p, x);
233 } else if (x < 0.0q) {
234 double_G_ibp (G, p, x);
235 } else {
236 double_G_cfrac_upper (G, p, x);
237 }
238 }
239
240 //! @brief iteration of the Romberg approximation of I_{x,y}^{mu,p}
241
242 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)
243 {
244 INT_T adr0_prev = ((n - 1) * n) / 2, adr0 = (n * (n + 1)) / 2;
245 DOUBLE_T sum = 0.0q;
246 for (INT_T j = 1; j <= pow2; j++) {
247 DOUBLE_T xx = x + ((y - x) * (2.0q * j - 1.0q)) / (2.0q * pow2);
248 sum += exp_double (-mu * xx + (p - 1.0q) * log_double (xx) - sigma);
249 }
250 R[adr0] = 0.5q * R[adr0_prev] + h * sum;
251 DOUBLE_T pow4 = 4.0q;
252 for (INT_T m = 1; m <= n; m++) {
253 R[adr0 + m] = (pow4 * R[adr0 + (m - 1)] - R[adr0_prev + (m - 1)]) / (pow4 - 1.0q);
254 pow4 *= 4.0q;
255 }
256 }
257
258 //! @ compute I_{x,y}^{mu,p} using a Romberg approximation.
259 // Compute rho and sigma so I_{x,y}^{mu,p} = rho * exp (sigma)
260
261 void double_romberg_estimate (DOUBLE_T * rho, DOUBLE_T * sigma, DOUBLE_T x, DOUBLE_T y, DOUBLE_T mu, DOUBLE_T p)
262 {
263 DOUBLE_T *R = (DOUBLE_T *) get_heap_space (((NITERMAX_ROMBERG + 1) * (NITERMAX_ROMBERG + 2)) / 2 * sizeof (DOUBLE_T));
264 ASSERT (R != NULL);
265 // Initialization (n=1)
266 *sigma = -mu * y + (p - 1.0q) * log_double (y);
267 R[0] = 0.5q * (y - x) * (exp_double (-mu * x + (p - 1.0q) * log_double (x) - (*sigma)) + 1.0q);
268 // Loop for n > 0
269 DOUBLE_T relneeded = A68_DOUBLE_EPS / TOL_ROMBERG;
270 INT_T adr0 = 0;
271 INT_T n = 1;
272 DOUBLE_T h = (y - x) / 2.0q; // n=1, h = (y-x)/2^n
273 DOUBLE_T pow2 = 1.0q; // n=1; pow2 = 2^(n-1)
274 if (NITERMAX_ROMBERG >= 1) {
275 DOUBLE_T relerr;
276 do {
277 double_romberg_iterations (R, *sigma, n, x, y, mu, p, h, pow2);
278 h /= 2.0q;
279 pow2 *= 2.0q;
280 adr0 = (n * (n + 1)) / 2;
281 relerr = fabs_double ((R[adr0 + n] - R[adr0 + n - 1]) / R[adr0 + n]);
282 n++;
283 } while (n <= NITERMAX_ROMBERG && relerr > relneeded);
284 }
285 // save Romberg estimate and free memory
286 *rho = R[adr0 + (n - 1)];
287 a68_free (R);
288 }
289
290 //! @brief compute generalized incomplete gamma function I_{x,y}^{mu,p}
291 // I_{x,y}^{mu,p} = integral from x to y of s^{p-1} * exp (-mu*s) ds
292 // This procedure computes (rho, sigma) described below.
293 // The approximated value of I_{x,y}^{mu,p} is I = rho * exp (sigma)
294 // mu is a real number non equal to zero
295 // (in general we take mu = 1 or -1 but any nonzero real number is allowed)
296 // x, y are two numbers with 0 <= x <= y <= +infinity,
297 // (the setting y=+infinity is allowed only when mu > 0)
298 // p is a real number > 0, p must be an integer when mu < 0.
299
300 void deltagammainc_double (DOUBLE_T * rho, DOUBLE_T * sigma, DOUBLE_T x, DOUBLE_T y, DOUBLE_T mu, DOUBLE_T p)
301 {
302 // Particular cases
303 if (isinf_double (x) && isinf_double (y)) {
304 *rho = 0.0q;
305 *sigma = a68_neginf_double ();
306 return;
307 } else if (x == y) {
308 *rho = 0.0q;
309 *sigma = a68_neginf_double ();
310 return;
311 }
312 if (x == 0.0q && isinf_double (y)) {
313 *rho = 1.0q;
314 (*sigma) = lgamma_double (p) - p * log_double (mu);
315 return;
316 }
317 // Initialization
318 DOUBLE_T mA, mB, mx, my, nA, nB;
319
320 double_G_func (&mx, p, mu * x);
321 DOUBLE_T nx = (isinf_double (x) ? a68_neginf_double () : -mu * x + p * log_double (x));
322 double_G_func (&my, p, mu * y);
323 DOUBLE_T ny = (isinf_double (y) ? a68_neginf_double () : -mu * y + p * log_double (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.0q) {
332 mA = my;
333 nA = ny;
334 mB = mx;
335 nB = nx;
336 } else {
337 if (p < double_plim (mu * x)) {
338 mA = mx;
339 nA = nx;
340 mB = my;
341 nB = ny;
342 } else if (p < double_plim (mu * y)) {
343 mA = 1.0q;
344 nA = lgamma_double (p) - p * log_double (mu);
345 nB = fmax (nx, ny);
346 mB = mx * exp_double (nx - nB) + my * exp_double (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 * exp_double (nB - nA);
356 *sigma = nA;
357 // If the difference involved a significant loss of precision, compute Romberg estimate.
358 if (!isinf_double (y) && ((*rho) / mA < TOL_DIFF)) {
359 double_romberg_estimate (rho, sigma, x, y, mu, p);
360 }
361 }
362
363 // A68G Driver routines
364
365 //! @brief PROC long gamma inc g = (LONG REAL p, x, y, mu) LONG REAL
366
367 void genie_gamma_inc_g_double (NODE_T * n)
368 {
369 A68_LONG_REAL x, y, mu, p;
370 POP_OBJECT (n, &mu, A68_LONG_REAL);
371 POP_OBJECT (n, &y, A68_LONG_REAL);
372 POP_OBJECT (n, &x, A68_LONG_REAL);
373 POP_OBJECT (n, &p, A68_LONG_REAL);
374 DOUBLE_T rho, sigma;
375 deltagammainc_double (&rho, &sigma, VALUE (&x).f, VALUE (&y).f, VALUE (&mu).f, VALUE (&p).f);
376 PUSH_VALUE (n, dble (rho * exp_double (sigma)), A68_LONG_REAL);
377 }
378
379 //! @brief PROC long gamma inc f = (LONG REAL p, x) LONG REAL
380
381 void genie_gamma_inc_f_double (NODE_T * n)
382 {
383 A68_LONG_REAL x, p;
384 POP_OBJECT (n, &x, A68_LONG_REAL);
385 POP_OBJECT (n, &p, A68_LONG_REAL);
386 DOUBLE_T rho, sigma;
387 deltagammainc_double (&rho, &sigma, VALUE (&x).f, a68_posinf_double (), 1.0q, VALUE (&p).f);
388 PUSH_VALUE (n, dble (rho * exp_double (sigma)), A68_LONG_REAL);
389 }
390
391 //! @brief PROC long gamma inc gf = (LONG REAL p, x) LONG REAL
392
393 void genie_gamma_inc_gf_double (NODE_T * q)
394 {
395 // if x <= p: G(p,x) = exp (x-p*ln (|x|)) * integral over [0,|x|] of s^{p-1} * exp (-sign (x)*s) ds
396 // otherwise: G(p,x) = exp (x-p*ln (x)) * integral over [x,inf] of s^{p-1} * exp (-s) ds
397 A68_LONG_REAL x, p;
398 POP_OBJECT (q, &x, A68_LONG_REAL);
399 POP_OBJECT (q, &p, A68_LONG_REAL);
400 DOUBLE_T G;
401 double_G_func (&G, VALUE (&p).f, VALUE (&x).f);
402 PUSH_VALUE (q, dble (G), A68_LONG_REAL);
403 }
404
405 //! @brief PROC long gamma inc = (LONG REAL p, x) LONG REAL
406
407 void genie_gamma_inc_h_double (NODE_T * n)
408 {
409 #if (A68_LEVEL >= 3) && defined (HAVE_GNU_MPFR)
410 genie_gamma_inc_double_mpfr (n);
411 #else
412 genie_gamma_inc_f_double (n);
413 #endif
414 }
415
416 #endif
