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-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 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 a68_max_real (REAL_T x, REAL_T y)
39 {
40 return (x > y ? x : y);
41 }
42
43 inline REAL_T a68_min_real (REAL_T x, REAL_T y)
44 {
45 return (x < y ? x : y);
46 }
47
48 inline REAL_T a68_sign_real (REAL_T x)
49 {
50 return (x == 0 ? 0 : (x > 0 ? 1 : -1));
51 }
52
53 inline REAL_T a68_int_real (REAL_T x)
54 {
55 return (x >= 0 ? (INT_T) x : -(INT_T) - x);
56 }
57
58 inline INT_T a68_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 == a68_int_real (n))
64
65 inline REAL_T a68_abs_real (REAL_T x)
66 {
67 return (x >= 0 ? x : -x);
68 }
69
70 REAL_T a68_fdiv_real (REAL_T x, REAL_T y)
71 {
72 // This is for generating +-INF.
73 return x / y;
74 }
75
76 REAL_T a68_nan_real (void)
77 {
78 return a68_fdiv_real (0.0, 0.0);
79 }
80
81 REAL_T a68_posinf_real (void)
82 {
83 return a68_fdiv_real (+1.0, 0.0);
84 }
85
86 REAL_T a68_neginf_double_real (void)
87 {
88 return a68_fdiv_real (-1.0, 0.0);
89 }
90
91 // REAL infinity
92
93 void genie_infinity_real (NODE_T * p)
94 {
95 PUSH_VALUE (p, a68_posinf_real (), A68_REAL);
96 }
97
98 // REAL minus infinity
99
100 void genie_minus_infinity_real (NODE_T * p)
101 {
102 PUSH_VALUE (p, a68_neginf_double_real (), A68_REAL);
103 }
104
105 int a68_finite_real (REAL_T x)
106 {
107 #if defined (HAVE_ISFINITE)
108 return isfinite (x);
109 #elif defined (HAVE_FINITE)
110 return finite (x);
111 #else
112 (void) x;
113 return A68_TRUE;
114 #endif
115 }
116
117 int a68_isnan_real (REAL_T x)
118 {
119 #if defined (HAVE_ISNAN)
120 return isnan (x);
121 #elif defined (HAVE_IEEE_COMPARISONS)
122 int status = (x != x);
123 return status;
124 #else
125 (void) x;
126 return A68_FALSE;
127 #endif
128 }
129
130 int a68_isinf_real (REAL_T x)
131 {
132 #if defined (HAVE_ISINF)
133 if (isinf (x)) {
134 return (x > 0) ? 1 : -1;
135 } else {
136 return 0;
137 }
138 #else
139 if (!a68_finite_real (x) && !a68_isnan_real (x)) {
140 return (x > 0 ? 1 : -1);
141 } else {
142 return 0;
143 }
144 #endif
145 }
146
147 // INT operators
148
149 INT_T a68_add_int (INT_T j, INT_T k)
150 {
151 if (j >= 0) {
152 A68_OVERFLOW (A68_MAX_INT - j < k);
153 } else {
154 A68_OVERFLOW (k < (-A68_MAX_INT) - j);
155 }
156 return j + k;
157 }
158
159 INT_T a68_sub_int (INT_T j, INT_T k)
160 {
161 return a68_add_int (j, -k);
162 }
163
164 INT_T a68_mul_int (INT_T j, INT_T k)
165 {
166 if (j == 0 || k == 0) {
167 return 0;
168 } else {
169 INT_T u = (j > 0 ? j : -j), v = (k > 0 ? k : -k);
170 A68_OVERFLOW (u > A68_MAX_INT / v);
171 return j * k;
172 }
173 }
174
175 INT_T a68_over_int (INT_T j, INT_T k)
176 {
177 A68_INVALID (k == 0);
178 return j / k;
179 }
180
181 INT_T a68_mod_int (INT_T j, INT_T k)
182 {
183 A68_INVALID (k == 0);
184 INT_T r = j % k;
185 return (r < 0 ? (k > 0 ? r + k : r - k) : r);
186 }
187
188 // OP ** = (INT, INT) INT
189
190 INT_T a68_m_up_n (INT_T m, INT_T n)
191 {
192 // Only positive n.
193 A68_INVALID (n < 0);
194 // Special cases.
195 if (m == 0 || m == 1) {
196 return m;
197 } else if (m == -1) {
198 return (EVEN (n) ? 1 : -1);
199 }
200 // General case with overflow check.
201 UNSIGNED_T bit = 1;
202 INT_T M = m, P = 1;
203 do {
204 if (n & bit) {
205 P = a68_mul_int (P, M);
206 }
207 bit <<= 1;
208 if (bit <= n) {
209 M = a68_mul_int (M, M);
210 }
211 } while (bit <= n);
212 return P;
213 }
214
215 // OP ** = (REAL, INT) REAL
216
217 REAL_T a68_x_up_n_real (REAL_T x, INT_T n)
218 {
219 // Only positive n.
220 if (n < 0) {
221 return 1 / a68_x_up_n_real (x, -n);
222 }
223 // Special cases.
224 if (x == 0 || x == 1) {
225 return x;
226 } else if (x == -1) {
227 return (EVEN (n) ? 1 : -1);
228 }
229 // General case.
230 UNSIGNED_T bit = 1;
231 REAL_T M = x, P = 1;
232 do {
233 if (n & bit) {
234 P *= M;
235 }
236 bit <<= 1;
237 if (bit <= n) {
238 M *= M;
239 }
240 } while (bit <= n);
241 A68_OVERFLOW (!a68_finite_real (P));
242 return P;
243 }
244
245 REAL_T a68_div_int (INT_T j, INT_T k)
246 {
247 A68_INVALID (k == 0);
248 return (REAL_T) j / (REAL_T) k;
249 }
250
251 // Sqrt (x^2 + y^2) that does not needlessly overflow.
252
253 REAL_T a68_hypot_real (REAL_T x, REAL_T y)
254 {
255 REAL_T xabs = ABS (x), yabs = ABS (y), min, max;
256 if (xabs < yabs) {
257 min = xabs;
258 max = yabs;
259 } else {
260 min = yabs;
261 max = xabs;
262 }
263 if (min == 0) {
264 return max;
265 } else {
266 REAL_T u = min / max;
267 return max * sqrt (1 + u * u);
268 }
269 }
270
271 //! @brief Compute Chebyshev series to requested accuracy.
272
273 REAL_T a68_chebyshev_real (REAL_T x, const REAL_T c[], REAL_T acc)
274 {
275 // Iteratively compute the recursive Chebyshev series.
276 // c[1..N] are coefficients, c[0] is N, and acc is relative accuracy.
277 acc *= A68_REAL_EPS;
278 if (acc < c[1]) {
279 diagnostic (A68_MATH_WARNING, A68 (f_entry), WARNING_MATH_ACCURACY, NULL);
280 }
281 INT_T N = a68_round (c[0]);
282 REAL_T err = 0, z = 2 * x, u = 0, v = 0, w = 0;
283 for (int i = 1; i <= N; i++) {
284 if (err > acc) {
285 w = v;
286 v = u;
287 u = z * v - w + c[i];
288 }
289 err += a68_abs_real (c[i]);
290 }
291 return 0.5 * (u - w);
292 }
293
294 // Compute ln (1 + x) accurately.
295 // Some C99 platforms implement this incorrectly.
296
297 REAL_T a68_ln1p_real (REAL_T x)
298 {
299 // Based on GNU GSL's gsl_sf_log_1plusx_e.
300 A68_INVALID (x <= -1);
301 if (a68_abs_real (x) < pow (A68_REAL_EPS, 1 / 6.0)) {
302 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;
303 const REAL_T t = c5 + x * (c6 + x * (c7 + x * (c8 + x * c9)));
304 return x * (1 + x * (c1 + x * (c2 + x * (c3 + x * (c4 + x * t)))));
305 } else if (a68_abs_real (x) < 0.5) {
306 REAL_T t = (8 * x + 1) / (x + 2) / 2;
307 return x * a68_chebyshev_real (t, c_ln1p, 0.1);
308 } else {
309 return ln (1 + x);
310 }
311 }
312
313 // Compute ln (x), if possible accurately when x ~ 1.
314
315 REAL_T a68_ln_real (REAL_T x)
316 {
317 A68_INVALID (x <= 0);
318 #if (A68_LEVEL >= 3)
319 if (a68_abs_real (x - 1) < 0.375) {
320 // Double precision x-1 mitigates cancellation error.
321 return a68_ln1p_real (DBLEQ (x) - 1.0q);
322 } else {
323 return ln (x);
324 }
325 #else
326 return ln (x);
327 #endif
328 }
329
330 // PROC (REAL) REAL exp
331
332 REAL_T a68_exp_real (REAL_T x)
333 {
334 return exp (x);
335 }
336
337 // OP ** = (REAL, REAL) REAL
338
339 REAL_T a68_x_up_y (REAL_T x, REAL_T y)
340 {
341 return a68_exp_real (y * a68_ln_real (x));
342 }
343
344 // PROC (REAL) REAL csc
345
346 REAL_T a68_csc_real (REAL_T x)
347 {
348 REAL_T z = sin (x);
349 A68_OVERFLOW (z == 0);
350 return 1 / z;
351 }
352
353 // PROC (REAL) REAL acsc
354
355 REAL_T a68_acsc_real (REAL_T x)
356 {
357 A68_OVERFLOW (x == 0);
358 return asin (1 / x);
359 }
360
361 // PROC (REAL) REAL sec
362
363 REAL_T a68_sec_real (REAL_T x)
364 {
365 REAL_T z = cos (x);
366 A68_OVERFLOW (z == 0);
367 return 1 / z;
368 }
369
370 // PROC (REAL) REAL asec
371
372 REAL_T a68_asec_real (REAL_T x)
373 {
374 A68_OVERFLOW (x == 0);
375 return acos (1 / x);
376 }
377
378 // PROC (REAL) REAL cot
379
380 REAL_T a68_cot_real (REAL_T x)
381 {
382 REAL_T z = sin (x);
383 A68_OVERFLOW (z == 0);
384 return cos (x) / z;
385 }
386
387 // PROC (REAL) REAL acot
388
389 REAL_T a68_acot_real (REAL_T x)
390 {
391 A68_OVERFLOW (x == 0);
392 return atan (1 / x);
393 }
394
395 // PROC atan2 (REAL, REAL) REAL
396
397 REAL_T a68_atan2_real (REAL_T x, REAL_T y)
398 {
399 if (x == 0) {
400 A68_INVALID (y == 0);
401 return (y > 0 ? CONST_PI_2 : -CONST_PI_2);
402 } else {
403 REAL_T z = atan (ABS (y / x));
404 if (x < 0) {
405 z = CONST_PI - z;
406 }
407 return (y >= 0 ? z : -z);
408 }
409 }
410
411 //! brief PROC (REAL) REAL cas
412
413 REAL_T a68_cas_real (REAL_T x)
414 {
415 // Hartley kernel, which Hartley named cas (cosine and sine).
416 return cos (x) + sin (x);
417 }
418
419 //! brief PROC (REAL) REAL sindg
420
421 REAL_T a68_sindg_real (REAL_T x)
422 {
423 return sin (x * CONST_PI_OVER_180);
424 }
425
426 //! brief PROC (REAL) REAL cosdg
427
428 REAL_T a68_cosdg_real (REAL_T x)
429 {
430 return cos (x * CONST_PI_OVER_180);
431 }
432
433 //! brief PROC (REAL) REAL tandg
434
435 REAL_T a68_tandg_real (REAL_T x)
436 {
437 return tan (x * CONST_PI_OVER_180);
438 }
439
440 //! brief PROC (REAL) REAL asindg
441
442 REAL_T a68_asindg_real (REAL_T x)
443 {
444 return asin (x) * CONST_180_OVER_PI;
445 }
446
447 //! brief PROC (REAL) REAL acosdg
448
449 REAL_T a68_acosdg_real (REAL_T x)
450 {
451 return acos (x) * CONST_180_OVER_PI;
452 }
453
454 //! brief PROC (REAL) REAL atandg
455
456 REAL_T a68_atandg_real (REAL_T x)
457 {
458 return atan (x) * CONST_180_OVER_PI;
459 }
460
461 // PROC (REAL) REAL cscd
462
463 REAL_T a68_cscdg_real (REAL_T x)
464 {
465 REAL_T z = a68_sindg_real (x);
466 A68_OVERFLOW (z == 0);
467 return 1 / z;
468 }
469
470 // PROC (REAL) REAL acscdg
471
472 REAL_T a68_acscdg_real (REAL_T x)
473 {
474 A68_OVERFLOW (x == 0);
475 return a68_asindg_real (1 / x);
476 }
477
478 // PROC (REAL) REAL secdg
479
480 REAL_T a68_secdg_real (REAL_T x)
481 {
482 REAL_T z = a68_cosdg_real (x);
483 A68_OVERFLOW (z == 0);
484 return 1 / z;
485 }
486
487 // PROC (REAL) REAL asecdg
488
489 REAL_T a68_asecdg_real (REAL_T x)
490 {
491 A68_OVERFLOW (x == 0);
492 return a68_acosdg_real (1 / x);
493 }
494
495 // PROC (REAL) REAL cotdg
496
497 REAL_T a68_cot_realdg_real (REAL_T x)
498 {
499 REAL_T z = a68_sindg_real (x);
500 A68_OVERFLOW (z == 0);
501 return a68_cosdg_real (x) / z;
502 }
503
504 // PROC (REAL) REAL acotdg
505
506 REAL_T a68_acotdg_real (REAL_T z)
507 {
508 A68_OVERFLOW (z == 0);
509 return a68_atandg_real (1 / z);
510 }
511
512 // @brief PROC (REAL) REAL sinpi
513
514 REAL_T a68_sinpi_real (REAL_T x)
515 {
516 x = fmod (x, 2);
517 if (x <= -1) {
518 x += 2;
519 } else if (x > 1) {
520 x -= 2;
521 }
522 // x in <-1, 1].
523 if (x == 0 || x == 1) {
524 return 0;
525 } else if (x == 0.5) {
526 return 1;
527 }
528 if (x == -0.5) {
529 return -1;
530 } else {
531 return sin (CONST_PI * x);
532 }
533 }
534
535 // @brief PROC (REAL) REAL cospi
536
537 REAL_T a68_cospi_real (REAL_T x)
538 {
539 x = fmod (fabs (x), 2);
540 // x in [0, 2>.
541 if (x == 0.5 || x == 1.5) {
542 return 0;
543 } else if (x == 0) {
544 return 1;
545 } else if (x == 1) {
546 return -1;
547 } else {
548 return cos (CONST_PI * x);
549 }
550 }
551
552 // @brief PROC (REAL) REAL tanpi
553
554 REAL_T a68_tanpi_real (REAL_T x)
555 {
556 x = fmod (x, 1);
557 if (x <= -0.5) {
558 x += 1;
559 } else if (x > 0.5) {
560 x -= 1;
561 }
562 // x in <-1/2, 1/2].
563 A68_OVERFLOW (x == 0.5);
564 if (x == -0.25) {
565 return -1;
566 } else if (x == 0) {
567 return 0;
568 } else if (x == 0.25) {
569 return 1;
570 } else {
571 return a68_sinpi_real (x) / a68_cospi_real (x);
572 }
573 }
574
575 // @brief PROC (REAL) REAL cotpi
576
577 REAL_T a68_cot_realpi (REAL_T x)
578 {
579 x = fmod (x, 1);
580 if (x <= -0.5) {
581 x += 1;
582 } else if (x > 0.5) {
583 x -= 1;
584 }
585 // x in <-1/2, 1/2].
586 A68_OVERFLOW (x == 0);
587 if (x == -0.25) {
588 return -1;
589 } else if (x == 0.25) {
590 return 1;
591 } else if (x == 0.5) {
592 return 0;
593 } else {
594 return a68_cospi_real (x) / a68_sinpi_real (x);
595 }
596 }
597
598 // @brief PROC (REAL) REAL asinh
599
600 REAL_T a68_asinh_real (REAL_T x)
601 {
602 REAL_T a = ABS (x), s = (x < 0 ? -1.0 : 1);
603 if (a > 1 / sqrt (A68_REAL_EPS)) {
604 return (s * (a68_ln_real (a) + a68_ln_real (2)));
605 } else if (a > 2) {
606 return (s * a68_ln_real (2 * a + 1 / (a + sqrt (a * a + 1))));
607 } else if (a > sqrt (A68_REAL_EPS)) {
608 REAL_T a2 = a * a;
609 return (s * a68_ln1p_real (a + a2 / (1 + sqrt (1 + a2))));
610 } else {
611 return (x);
612 }
613 }
614
615 // @brief PROC (REAL) REAL acosh
616
617 REAL_T a68_acosh_real (REAL_T x)
618 {
619 if (x > 1 / sqrt (A68_REAL_EPS)) {
620 return (a68_ln_real (x) + a68_ln_real (2));
621 } else if (x > 2) {
622 return (a68_ln_real (2 * x - 1 / (sqrt (x * x - 1) + x)));
623 } else if (x > 1) {
624 REAL_T t = x - 1;
625 return (a68_ln1p_real (t + sqrt (2 * t + t * t)));
626 } else if (x == 1) {
627 return (0);
628 } else {
629 A68_INVALID (A68_TRUE);
630 }
631 }
632
633 // @brief PROC (REAL) REAL atanh
634
635 REAL_T a68_atanh_real (REAL_T x)
636 {
637 REAL_T a = ABS (x);
638 A68_INVALID (a >= 1);
639 REAL_T s = (x < 0 ? -1 : 1);
640 if (a >= 0.5) {
641 return (s * 0.5 * a68_ln1p_real (2 * a / (1 - a)));
642 } else if (a > A68_REAL_EPS) {
643 return (s * 0.5 * a68_ln1p_real (2 * a + 2 * a * a / (1 - a)));
644 } else {
645 return (x);
646 }
647 }
648
649 //! @brief Inverse complementary error function.
650
651 REAL_T a68_inverfc_real (REAL_T y)
652 {
653 A68_INVALID (y < 0 || y > 2);
654 if (y == 0) {
655 return A68_REAL_MAX;
656 } else if (y == 1) {
657 return 0;
658 } else if (y == 2) {
659 return -A68_REAL_MAX;
660 } else {
661 // Next is based on code that originally contained following statement:
662 // Copyright (c) 1996 Takuya Ooura. You may use, copy, modify this
663 // code for any purpose and without fee.
664 REAL_T s, t, u, v, x, z;
665 z = (y <= 1 ? y : 2 - y);
666 v = c_inverfc[0] - a68_ln_real (z);
667 u = sqrt (v);
668 s = (a68_ln_real (u) + c_inverfc[1]) / v;
669 t = 1 / (u + c_inverfc[2]);
670 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;
671 t = c_inverfc[9] / (x + c_inverfc[9]);
672 u = t - 0.5;
673 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];
674 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 * a68_exp_real (x * x - c_inverfc[33]);
675 x += s * (x * s + 1);
676 return (y <= 1 ? x : -x);
677 }
678 }
679
680 //! @brief Inverse error function.
681
682 REAL_T a68_inverf_real (REAL_T y)
683 {
684 return a68_inverfc_real (1 - y);
685 }
686
687 //! @brief PROC (REAL, REAL) REAL ln beta
688
689 REAL_T a68_ln_beta_real (REAL_T a, REAL_T b)
690 {
691 return lgamma (a) + lgamma (b) - lgamma (a + b);
692 }
693
694 //! @brief PROC (REAL, REAL) REAL beta
695
696 REAL_T a68_beta_real (REAL_T a, REAL_T b)
697 {
698 return a68_exp_real (a68_ln_beta_real (a, b));
699 }
700
701 //! brief PROC (INT) REAL fact
702
703 REAL_T a68_fact_real (INT_T n)
704 {
705 A68_INVALID (n < 0 || n > A68_MAX_FAC);
706 return factable[n];
707 }
708
709 //! brief PROC (INT) REAL ln fact
710
711 REAL_T a68_ln_fact_real (INT_T n)
712 {
713 A68_INVALID (n < 0);
714 if (n <= A68_MAX_FAC) {
715 return ln_factable[n];
716 } else {
717 return lgamma (n + 1);
718 }
719 }
720
721 //! @brief PROC choose = (INT n, m) REAL
722
723 REAL_T a68_choose_real (INT_T n, INT_T m)
724 {
725 A68_INVALID (n < m);
726 return factable[n] / (factable[m] * factable[n - m]);
727 }
728
729 //! @brief PROC ln choose = (INT n, m) REAL
730
731 REAL_T a68_ln_choose_real (INT_T n, INT_T m)
732 {
733 A68_INVALID (n < m);
734 return a68_ln_fact_real (n) - (a68_ln_fact_real (m) + a68_ln_fact_real (n - m));
735 }
736
737 REAL_T a68_beta_inc_real (REAL_T s, REAL_T t, REAL_T x)
738 {
739 // Incomplete beta function I{x}(s, t).
740 // Continued fraction, see dlmf.nist.gov/8.17; Lentz's algorithm.
741 if (x < 0 || x > 1) {
742 errno = ERANGE;
743 return -1;
744 } else {
745 const INT_T lim = 16 * sizeof (REAL_T);
746 BOOL_T cont = A68_TRUE;
747 // Rapid convergence when x <= (s+1)/(s+t+2) or else recursion.
748 if (x > (s + 1) / (s + t + 2)) {
749 // B{x}(s, t) = 1 - B{1-x}(t, s)
750 return 1 - a68_beta_inc_real (s, t, 1 - x);
751 }
752 // Lentz's algorithm for continued fraction.
753 REAL_T W = 1, F = 1, c = 1, d = 0;
754 INT_T N, m;
755 for (N = 0, m = 0; cont && N < lim; N++) {
756 REAL_T T;
757 if (N == 0) {
758 T = 1;
759 } else if (N % 2 == 0) {
760 // d{2m} := x m(t-m)/((s+2m-1)(s+2m))
761 T = x * m * (t - m) / (s + 2 * m - 1) / (s + 2 * m);
762 } else {
763 // d{2m+1} := -x (s+m)(s+t+m)/((s+2m+1)(s+2m))
764 T = -x * (s + m) * (s + t + m) / (s + 2 * m + 1) / (s + 2 * m);
765 m++;
766 }
767 d = 1 / (T * d + 1);
768 c = T / c + 1;
769 F *= c * d;
770 if (F == W) {
771 cont = A68_FALSE;
772 } else {
773 W = F;
774 }
775 }
776 // I{x}(s,t)=x^s(1-x)^t / s / B(s,t) F
777 REAL_T beta = a68_exp_real (lgamma (s) + lgamma (t) - lgamma (s + t));
778 return a68_x_up_y (x, s) * a68_x_up_y (1 - x, t) / s / beta * (F - 1);
779 }
780 }
