single-math.c
1 //! @file single-math.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-2025 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 math routines supplementing libc.
25
26 // References:
27 // M. Abramowitz and I. Stegun, Handbook of Mathematical Functions,
28 // Dover Publications, New York [1970]
29 // https://en.wikipedia.org/wiki/Abramowitz_and_Stegun
30
31 #include "a68g.h"
32 #include "a68g-genie.h"
33 #include "a68g-prelude.h"
34 #include "a68g-double.h"
35 #include "a68g-numbers.h"
36 #include "a68g-math.h"
37
38 inline REAL_T a68g_max_real (REAL_T x, REAL_T y)
39 {
40 return (x > y ? x : y);
41 }
42
43 inline REAL_T a68g_min_real (REAL_T x, REAL_T y)
44 {
45 return (x < y ? x : y);
46 }
47
48 inline REAL_T a68g_sign_real (REAL_T x)
49 {
50 return (x == 0 ? 0 : (x > 0 ? 1 : -1));
51 }
52
53 inline REAL_T a68g_int_real (REAL_T x)
54 {
55 return (x >= 0 ? (INT_T) x : -(INT_T) - x);
56 }
57
58 inline INT_T a68g_round (REAL_T x)
59 {
60 return (INT_T) (x >= 0 ? x + 0.5 : x - 0.5);
61 }
62
63 #define IS_INTEGER(n) (n == a68g_int_real (n))
64
65 inline REAL_T a68g_abs_real (REAL_T x)
66 {
67 return (x >= 0 ? x : -x);
68 }
69
70 REAL_T a68g_fdiv_real (REAL_T x, REAL_T y)
71 {
72 // This is for generating +-INF.
73 return x / y;
74 }
75
76 REAL_T a68g_nan_real (void)
77 {
78 return a68g_fdiv_real (0.0, 0.0);
79 }
80
81 REAL_T a68g_posinf_real (void)
82 {
83 return a68g_fdiv_real (+1.0, 0.0);
84 }
85
86 REAL_T a68g_mininf_real (void)
87 {
88 return a68g_fdiv_real (-1.0, 0.0);
89 }
90
91 // Range check.
92
93 BOOL_T a68g_finite_real (REAL_T x)
94 {
95 #if defined (HAVE_ISFINITE)
96 return isfinite (x);
97 #elif defined (HAVE_FINITE)
98 return finite (x);
99 #else
100 (void) x;
101 return A68G_TRUE;
102 #endif
103 }
104
105 // Range check.
106
107 BOOL_T a68g_isnan_real (REAL_T x)
108 {
109 #if defined (HAVE_ISNAN)
110 return isnan (x) == FP_NAN;
111 #elif defined (HAVE_IEEE_COMPARISONS)
112 int status = (x != x);
113 return status;
114 #else
115 (void) x;
116 return A68G_FALSE;
117 #endif
118 }
119
120 // Range check.
121
122 BOOL_T a68g_isinf_real (REAL_T x)
123 {
124 #if defined (HAVE_ISINF)
125 return isinf (x);
126 #else
127 if (!a68g_finite_real (x) && !a68g_isnan_real (x)) {
128 return A68G_TRUE
129 } else {
130 return A68G_FALSE
131 }
132 #endif
133 }
134
135 // INT operators
136
137 INT_T a68g_add_int (INT_T j, INT_T k)
138 {
139 if (j >= 0) {
140 A68G_OVERFLOW (A68G_MAX_INT - j < k);
141 } else {
142 A68G_OVERFLOW (k < (-A68G_MAX_INT) - j);
143 }
144 return j + k;
145 }
146
147 INT_T a68g_sub_int (INT_T j, INT_T k)
148 {
149 return a68g_add_int (j, -k);
150 }
151
152 INT_T a68g_mul_int (INT_T j, INT_T k)
153 {
154 if (j == 0 || k == 0) {
155 return 0;
156 } else {
157 INT_T u = (j > 0 ? j : -j), v = (k > 0 ? k : -k);
158 A68G_OVERFLOW (u > A68G_MAX_INT / v);
159 return j * k;
160 }
161 }
162
163 INT_T a68g_over_int (INT_T j, INT_T k)
164 {
165 A68G_INVALID (k == 0);
166 return j / k;
167 }
168
169 INT_T a68g_mod_int (INT_T j, INT_T k)
170 {
171 A68G_INVALID (k == 0);
172 INT_T r = j % k;
173 return (r < 0 ? (k > 0 ? r + k : r - k) : r);
174 }
175
176 // OP ** = (INT, INT) INT
177
178 INT_T a68g_m_up_n (INT_T m, INT_T n)
179 {
180 // Only positive n.
181 A68G_INVALID (n < 0);
182 // Special cases.
183 if (m == 0 || m == 1) {
184 return m;
185 } else if (m == -1) {
186 return (EVEN (n) ? 1 : -1);
187 }
188 // General case with overflow check.
189 UNSIGNED_T bit = 1;
190 INT_T M = m, P = 1;
191 do {
192 if (n & bit) {
193 P = a68g_mul_int (P, M);
194 }
195 bit <<= 1;
196 if (bit <= n) {
197 M = a68g_mul_int (M, M);
198 }
199 } while (bit <= n);
200 return P;
201 }
202
203 // OP ** = (REAL, INT) REAL
204
205 REAL_T a68g_x_up_n_real (REAL_T x, INT_T n)
206 {
207 // Only positive n.
208 if (n < 0) {
209 return 1 / a68g_x_up_n_real (x, -n);
210 }
211 // Special cases.
212 if (x == 0 || x == 1) {
213 return x;
214 } else if (x == -1) {
215 return (EVEN (n) ? 1 : -1);
216 }
217 // General case.
218 UNSIGNED_T bit = 1;
219 REAL_T M = x, P = 1;
220 do {
221 if (n & bit) {
222 P *= M;
223 }
224 bit <<= 1;
225 if (bit <= n) {
226 M *= M;
227 }
228 } while (bit <= n);
229 A68G_OVERFLOW (!a68g_finite_real (P));
230 return P;
231 }
232
233 REAL_T a68g_div_int (INT_T j, INT_T k)
234 {
235 A68G_INVALID (k == 0);
236 return (REAL_T) j / (REAL_T) k;
237 }
238
239 // Sqrt (x^2 + y^2) that does not needlessly overflow.
240
241 REAL_T a68g_hypot_real (REAL_T x, REAL_T y)
242 {
243 REAL_T xabs = ABS (x), yabs = ABS (y), min, max;
244 if (xabs < yabs) {
245 min = xabs;
246 max = yabs;
247 } else {
248 min = yabs;
249 max = xabs;
250 }
251 if (min == 0) {
252 return max;
253 } else {
254 REAL_T u = min / max;
255 return max * sqrt (1 + u * u);
256 }
257 }
258
259 //! @brief Compute Chebyshev series to requested accuracy.
260
261 REAL_T a68g_chebyshev_real (REAL_T x, const REAL_T c[], REAL_T acc)
262 {
263 // Iteratively compute the recursive Chebyshev series.
264 // c[1..N] are coefficients, c[0] is N, and acc is relative accuracy.
265 acc *= A68G_REAL_EPS;
266 if (acc < c[1]) {
267 diagnostic (A68G_MATH_WARNING, A68G (f_entry), WARNING_MATH_ACCURACY, NULL);
268 }
269 INT_T N = a68g_round (c[0]);
270 REAL_T err = 0, z = 2 * x, u = 0, v = 0, w = 0;
271 for (int i = 1; i <= N; i++) {
272 if (err > acc) {
273 w = v;
274 v = u;
275 u = z * v - w + c[i];
276 }
277 err += a68g_abs_real (c[i]);
278 }
279 return 0.5 * (u - w);
280 }
281
282 // Compute ln (1 + x) accurately.
283 // Some C99 platforms implement this incorrectly.
284
285 REAL_T a68g_ln1p_real (REAL_T x)
286 {
287 // Based on GNU GSL's gsl_sf_log_1plusx_e.
288 A68G_INVALID (x <= -1);
289 if (a68g_abs_real (x) < pow (A68G_REAL_EPS, 1 / 6.0)) {
290 const REAL_T c1 = -0.5, c2 = 1 / 3.0, c3 = -1 / 4.0, c4 = 1 / 5.0, c5 = -1 / 6.0, c6 = 1 / 7.0, c7 = -1 / 8.0, c8 = 1 / 9.0, c9 = -1 / 10.0;
291 const REAL_T t = c5 + x * (c6 + x * (c7 + x * (c8 + x * c9)));
292 return x * (1 + x * (c1 + x * (c2 + x * (c3 + x * (c4 + x * t)))));
293 } else if (a68g_abs_real (x) < 0.5) {
294 REAL_T t = (8 * x + 1) / (x + 2) / 2;
295 return x * a68g_chebyshev_real (t, c_ln1p, 0.1);
296 } else {
297 return ln (1 + x);
298 }
299 }
300
301 // Compute ln (x), if possible accurately when x ~ 1.
302
303 REAL_T a68g_ln_real (REAL_T x)
304 {
305 A68G_INVALID (x <= 0);
306 #if (A68G_LEVEL >= 3)
307 if (a68g_abs_real (x - 1) < 0.375) {
308 // Double precision x-1 mitigates cancellation error.
309 return a68g_ln1p_real (DBLEQ (x) - 1.0q);
310 } else {
311 return ln (x);
312 }
313 #else
314 return ln (x);
315 #endif
316 }
317
318 // PROC (REAL) REAL exp
319
320 REAL_T a68g_exp_real (REAL_T x)
321 {
322 return exp (x);
323 }
324
325 // OP ** = (REAL, REAL) REAL
326
327 REAL_T a68g_x_up_y (REAL_T x, REAL_T y)
328 {
329 return a68g_exp_real (y * a68g_ln_real (x));
330 }
331
332 // PROC (REAL) REAL csc
333
334 REAL_T a68g_csc_real (REAL_T x)
335 {
336 REAL_T z = sin (x);
337 A68G_OVERFLOW (z == 0);
338 return 1 / z;
339 }
340
341 // PROC (REAL) REAL acsc
342
343 REAL_T a68g_acsc_real (REAL_T x)
344 {
345 A68G_OVERFLOW (x == 0);
346 return asin (1 / x);
347 }
348
349 // PROC (REAL) REAL sec
350
351 REAL_T a68g_sec_real (REAL_T x)
352 {
353 REAL_T z = cos (x);
354 A68G_OVERFLOW (z == 0);
355 return 1 / z;
356 }
357
358 // PROC (REAL) REAL asec
359
360 REAL_T a68g_asec_real (REAL_T x)
361 {
362 A68G_OVERFLOW (x == 0);
363 return acos (1 / x);
364 }
365
366 // PROC (REAL) REAL cot
367
368 REAL_T a68g_cot_real (REAL_T x)
369 {
370 REAL_T z = sin (x);
371 A68G_OVERFLOW (z == 0);
372 return cos (x) / z;
373 }
374
375 // PROC (REAL) REAL acot
376
377 REAL_T a68g_acot_real (REAL_T x)
378 {
379 A68G_OVERFLOW (x == 0);
380 return atan (1 / x);
381 }
382
383 // PROC atan2 (REAL, REAL) REAL
384
385 REAL_T a68g_atan2_real (REAL_T x, REAL_T y)
386 {
387 if (x == 0) {
388 A68G_INVALID (y == 0);
389 return (y > 0 ? CONST_PI_2 : -CONST_PI_2);
390 } else {
391 REAL_T z = atan (ABS (y / x));
392 if (x < 0) {
393 z = CONST_PI - z;
394 }
395 return (y >= 0 ? z : -z);
396 }
397 }
398
399 //! brief PROC (REAL) REAL cas
400
401 REAL_T a68g_cas_real (REAL_T x)
402 {
403 // Hartley kernel, which Hartley named cas (cosine and sine).
404 return cos (x) + sin (x);
405 }
406
407 //! brief PROC (REAL) REAL sindg
408
409 REAL_T a68g_sindg_real (REAL_T x)
410 {
411 return sin (x * CONST_PI_OVER_180);
412 }
413
414 //! brief PROC (REAL) REAL cosdg
415
416 REAL_T a68g_cosdg_real (REAL_T x)
417 {
418 return cos (x * CONST_PI_OVER_180);
419 }
420
421 //! brief PROC (REAL) REAL tandg
422
423 REAL_T a68g_tandg_real (REAL_T x)
424 {
425 return tan (x * CONST_PI_OVER_180);
426 }
427
428 //! brief PROC (REAL) REAL asindg
429
430 REAL_T a68g_asindg_real (REAL_T x)
431 {
432 return asin (x) * CONST_180_OVER_PI;
433 }
434
435 //! brief PROC (REAL) REAL acosdg
436
437 REAL_T a68g_acosdg_real (REAL_T x)
438 {
439 return acos (x) * CONST_180_OVER_PI;
440 }
441
442 //! brief PROC (REAL) REAL atandg
443
444 REAL_T a68g_atandg_real (REAL_T x)
445 {
446 return atan (x) * CONST_180_OVER_PI;
447 }
448
449 // PROC (REAL) REAL cscd
450
451 REAL_T a68g_cscdg_real (REAL_T x)
452 {
453 REAL_T z = a68g_sindg_real (x);
454 A68G_OVERFLOW (z == 0);
455 return 1 / z;
456 }
457
458 // PROC (REAL) REAL acscdg
459
460 REAL_T a68g_acscdg_real (REAL_T x)
461 {
462 A68G_OVERFLOW (x == 0);
463 return a68g_asindg_real (1 / x);
464 }
465
466 // PROC (REAL) REAL secdg
467
468 REAL_T a68g_secdg_real (REAL_T x)
469 {
470 REAL_T z = a68g_cosdg_real (x);
471 A68G_OVERFLOW (z == 0);
472 return 1 / z;
473 }
474
475 // PROC (REAL) REAL asecdg
476
477 REAL_T a68g_asecdg_real (REAL_T x)
478 {
479 A68G_OVERFLOW (x == 0);
480 return a68g_acosdg_real (1 / x);
481 }
482
483 // PROC (REAL) REAL cotdg
484
485 REAL_T a68g_cot_realdg_real (REAL_T x)
486 {
487 REAL_T z = a68g_sindg_real (x);
488 A68G_OVERFLOW (z == 0);
489 return a68g_cosdg_real (x) / z;
490 }
491
492 // PROC (REAL) REAL acotdg
493
494 REAL_T a68g_acotdg_real (REAL_T z)
495 {
496 A68G_OVERFLOW (z == 0);
497 return a68g_atandg_real (1 / z);
498 }
499
500 // @brief PROC (REAL) REAL sinpi
501
502 REAL_T a68g_sinpi_real (REAL_T x)
503 {
504 x = fmod (x, 2);
505 if (x <= -1) {
506 x += 2;
507 } else if (x > 1) {
508 x -= 2;
509 }
510 // x in <-1, 1].
511 if (x == 0 || x == 1) {
512 return 0;
513 } else if (x == 0.5) {
514 return 1;
515 }
516 if (x == -0.5) {
517 return -1;
518 } else {
519 return sin (CONST_PI * x);
520 }
521 }
522
523 // @brief PROC (REAL) REAL cospi
524
525 REAL_T a68g_cospi_real (REAL_T x)
526 {
527 x = fmod (fabs (x), 2);
528 // x in [0, 2>.
529 if (x == 0.5 || x == 1.5) {
530 return 0;
531 } else if (x == 0) {
532 return 1;
533 } else if (x == 1) {
534 return -1;
535 } else {
536 return cos (CONST_PI * x);
537 }
538 }
539
540 // @brief PROC (REAL) REAL tanpi
541
542 REAL_T a68g_tanpi_real (REAL_T x)
543 {
544 x = fmod (x, 1);
545 if (x <= -0.5) {
546 x += 1;
547 } else if (x > 0.5) {
548 x -= 1;
549 }
550 // x in <-1/2, 1/2].
551 A68G_OVERFLOW (x == 0.5);
552 if (x == -0.25) {
553 return -1;
554 } else if (x == 0) {
555 return 0;
556 } else if (x == 0.25) {
557 return 1;
558 } else {
559 return a68g_sinpi_real (x) / a68g_cospi_real (x);
560 }
561 }
562
563 // @brief PROC (REAL) REAL cotpi
564
565 REAL_T a68g_cot_realpi (REAL_T x)
566 {
567 x = fmod (x, 1);
568 if (x <= -0.5) {
569 x += 1;
570 } else if (x > 0.5) {
571 x -= 1;
572 }
573 // x in <-1/2, 1/2].
574 A68G_OVERFLOW (x == 0);
575 if (x == -0.25) {
576 return -1;
577 } else if (x == 0.25) {
578 return 1;
579 } else if (x == 0.5) {
580 return 0;
581 } else {
582 return a68g_cospi_real (x) / a68g_sinpi_real (x);
583 }
584 }
585
586 // @brief PROC (REAL) REAL asinh
587
588 REAL_T a68g_asinh_real (REAL_T x)
589 {
590 REAL_T a = ABS (x), s = (x < 0 ? -1.0 : 1);
591 if (a > 1 / sqrt (A68G_REAL_EPS)) {
592 return (s * (a68g_ln_real (a) + a68g_ln_real (2)));
593 } else if (a > 2) {
594 return (s * a68g_ln_real (2 * a + 1 / (a + sqrt (a * a + 1))));
595 } else if (a > sqrt (A68G_REAL_EPS)) {
596 REAL_T a2 = a * a;
597 return (s * a68g_ln1p_real (a + a2 / (1 + sqrt (1 + a2))));
598 } else {
599 return (x);
600 }
601 }
602
603 // @brief PROC (REAL) REAL acosh
604
605 REAL_T a68g_acosh_real (REAL_T x)
606 {
607 if (x > 1 / sqrt (A68G_REAL_EPS)) {
608 return (a68g_ln_real (x) + a68g_ln_real (2));
609 } else if (x > 2) {
610 return (a68g_ln_real (2 * x - 1 / (sqrt (x * x - 1) + x)));
611 } else if (x > 1) {
612 REAL_T t = x - 1;
613 return (a68g_ln1p_real (t + sqrt (2 * t + t * t)));
614 } else if (x == 1) {
615 return (0);
616 } else {
617 A68G_INVALID (A68G_TRUE);
618 }
619 }
620
621 // @brief PROC (REAL) REAL atanh
622
623 REAL_T a68g_atanh_real (REAL_T x)
624 {
625 REAL_T a = ABS (x);
626 A68G_INVALID (a >= 1);
627 REAL_T s = (x < 0 ? -1 : 1);
628 if (a >= 0.5) {
629 return (s * 0.5 * a68g_ln1p_real (2 * a / (1 - a)));
630 } else if (a > A68G_REAL_EPS) {
631 return (s * 0.5 * a68g_ln1p_real (2 * a + 2 * a * a / (1 - a)));
632 } else {
633 return (x);
634 }
635 }
636
637 //! @brief Inverse complementary error function.
638
639 REAL_T a68g_inverfc_real (REAL_T y)
640 {
641 A68G_INVALID (y < 0 || y > 2);
642 if (y == 0) {
643 return A68G_REAL_MAX;
644 } else if (y == 1) {
645 return 0;
646 } else if (y == 2) {
647 return -A68G_REAL_MAX;
648 } else {
649 // Next is based on code that originally contained following statement:
650 // Copyright (c) 1996 Takuya Ooura. You may use, copy, modify this
651 // code for any purpose and without fee.
652 REAL_T s, t, u, v, x, z;
653 z = (y <= 1 ? y : 2 - y);
654 v = c_inverfc[0] - a68g_ln_real (z);
655 u = sqrt (v);
656 s = (a68g_ln_real (u) + c_inverfc[1]) / v;
657 t = 1 / (u + c_inverfc[2]);
658 x = u * (1 - s * (s * c_inverfc[3] + 0.5)) - ((((c_inverfc[4] * t + c_inverfc[5]) * t + c_inverfc[6]) * t + c_inverfc[7]) * t + c_inverfc[8]) * t;
659 t = c_inverfc[9] / (x + c_inverfc[9]);
660 u = t - 0.5;
661 s = (((((((((c_inverfc[10] * u + c_inverfc[11]) * u - c_inverfc[12]) * u - c_inverfc[13]) * u + c_inverfc[14]) * u + c_inverfc[15]) * u - c_inverfc[16]) * u - c_inverfc[17]) * u + c_inverfc[18]) * u + c_inverfc[19]) * u + c_inverfc[20];
662 s = ((((((((((((s * u - c_inverfc[21]) * u - c_inverfc[22]) * u + c_inverfc[23]) * u + c_inverfc[24]) * u + c_inverfc[25]) * u + c_inverfc[26]) * u + c_inverfc[27]) * u + c_inverfc[28]) * u + c_inverfc[29]) * u + c_inverfc[30]) * u + c_inverfc[31]) * u + c_inverfc[32]) * t - z * a68g_exp_real (x * x - c_inverfc[33]);
663 x += s * (x * s + 1);
664 return (y <= 1 ? x : -x);
665 }
666 }
667
668 //! @brief Inverse error function.
669
670 REAL_T a68g_inverf_real (REAL_T y)
671 {
672 return a68g_inverfc_real (1 - y);
673 }
674
675 //! @brief PROC (REAL, REAL) REAL ln beta
676
677 REAL_T a68g_ln_beta_real (REAL_T a, REAL_T b)
678 {
679 return lgamma (a) + lgamma (b) - lgamma (a + b);
680 }
681
682 //! @brief PROC (REAL, REAL) REAL beta
683
684 REAL_T a68g_beta_real (REAL_T a, REAL_T b)
685 {
686 return a68g_exp_real (a68g_ln_beta_real (a, b));
687 }
688
689 //! brief PROC (INT) REAL fact
690
691 REAL_T a68g_fact_real (INT_T n)
692 {
693 A68G_INVALID (n < 0 || n > A68G_MAX_FAC);
694 return factable[n];
695 }
696
697 //! brief PROC (INT) REAL ln fact
698
699 REAL_T a68g_ln_fact_real (INT_T n)
700 {
701 A68G_INVALID (n < 0);
702 if (n <= A68G_MAX_FAC) {
703 return ln_factable[n];
704 } else {
705 return lgamma (n + 1);
706 }
707 }
708
709 //! @brief PROC choose = (INT n, m) REAL
710
711 REAL_T a68g_choose_real (INT_T n, INT_T m)
712 {
713 A68G_INVALID (n < m);
714 return factable[n] / (factable[m] * factable[n - m]);
715 }
716
717 //! @brief PROC ln choose = (INT n, m) REAL
718
719 REAL_T a68g_ln_choose_real (INT_T n, INT_T m)
720 {
721 A68G_INVALID (n < m);
722 return a68g_ln_fact_real (n) - (a68g_ln_fact_real (m) + a68g_ln_fact_real (n - m));
723 }
724
725 REAL_T a68g_beta_inc_real (REAL_T s, REAL_T t, REAL_T x)
726 {
727 // Incomplete beta function I{x}(s, t).
728 // Continued fraction, see dlmf.nist.gov/8.17; Lentz's algorithm.
729 if (x < 0 || x > 1) {
730 errno = ERANGE;
731 return -1;
732 } else {
733 const INT_T lim = 16 * sizeof (REAL_T);
734 BOOL_T cont = A68G_TRUE;
735 // Rapid convergence when x <= (s+1)/(s+t+2) or else recursion.
736 if (x > (s + 1) / (s + t + 2)) {
737 // B{x}(s, t) = 1 - B{1-x}(t, s)
738 return 1 - a68g_beta_inc_real (s, t, 1 - x);
739 }
740 // Lentz's algorithm for continued fraction.
741 REAL_T W = 1, F = 1, c = 1, d = 0;
742 INT_T N, m;
743 for (N = 0, m = 0; cont && N < lim; N++) {
744 REAL_T T;
745 if (N == 0) {
746 T = 1;
747 } else if (N % 2 == 0) {
748 // d{2m} := x m(t-m)/((s+2m-1)(s+2m))
749 T = x * m * (t - m) / (s + 2 * m - 1) / (s + 2 * m);
750 } else {
751 // d{2m+1} := -x (s+m)(s+t+m)/((s+2m+1)(s+2m))
752 T = -x * (s + m) * (s + t + m) / (s + 2 * m + 1) / (s + 2 * m);
753 m++;
754 }
755 d = 1 / (T * d + 1);
756 c = T / c + 1;
757 F *= c * d;
758 if (F == W) {
759 cont = A68G_FALSE;
760 } else {
761 W = F;
762 }
763 }
764 // I{x}(s,t)=x^s(1-x)^t / s / B(s,t) F
765 REAL_T beta = a68g_exp_real (lgamma (s) + lgamma (t) - lgamma (s + t));
766 return a68g_x_up_y (x, s) * a68g_x_up_y (1 - x, t) / s / beta * (F - 1);
767 }
768 }
|
© 2002-2025 J.M. van der Veer (jmvdveer@xs4all.nl)
|
