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) != 0;
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 && n == 0) {
184 return 1;
185 } else if (m == 0 || m == 1) {
186 return m;
187 } else if (m == -1) {
188 return (EVEN (n) ? 1 : -1);
189 }
190 // General case with overflow check.
191 UNSIGNED_T bit = 1;
192 INT_T M = m, P = 1;
193 do {
194 if (n & bit) {
195 P = a68g_mul_int (P, M);
196 }
197 bit <<= 1;
198 if (bit <= n) {
199 M = a68g_mul_int (M, M);
200 }
201 } while (bit <= n);
202 return P;
203 }
204
205 // OP ** = (REAL, INT) REAL
206
207 REAL_T a68g_x_up_n_real (REAL_T x, INT_T n)
208 {
209 // Only positive n.
210 if (n < 0) {
211 return 1 / a68g_x_up_n_real (x, -n);
212 }
213 // Special cases.
214 if (x == 0 && n == 0) {
215 return 1;
216 } else if (x == 0 || x == 1) {
217 return x;
218 } else if (x == -1) {
219 return (EVEN (n) ? 1 : -1);
220 }
221 // General case.
222 UNSIGNED_T bit = 1;
223 REAL_T M = x, P = 1;
224 do {
225 if (n & bit) {
226 P *= M;
227 }
228 bit <<= 1;
229 if (bit <= n) {
230 M *= M;
231 }
232 } while (bit <= n);
233 A68G_OVERFLOW (!a68g_finite_real (P));
234 return P;
235 }
236
237 REAL_T a68g_div_int (INT_T j, INT_T k)
238 {
239 A68G_INVALID (k == 0);
240 return (REAL_T) j / (REAL_T) k;
241 }
242
243 // Sqrt (x^2 + y^2) that does not needlessly overflow.
244
245 REAL_T a68g_hypot_real (REAL_T x, REAL_T y)
246 {
247 REAL_T xabs = ABS (x), yabs = ABS (y), min, max;
248 if (xabs < yabs) {
249 min = xabs;
250 max = yabs;
251 } else {
252 min = yabs;
253 max = xabs;
254 }
255 if (min == 0) {
256 return max;
257 } else {
258 REAL_T u = min / max;
259 return max * sqrt (1 + u * u);
260 }
261 }
262
263 //! @brief Compute Chebyshev series to requested accuracy.
264
265 REAL_T a68g_chebyshev_real (REAL_T x, const REAL_T c[], REAL_T acc)
266 {
267 // Iteratively compute the recursive Chebyshev series.
268 // c[1..N] are coefficients, c[0] is N, and acc is relative accuracy.
269 acc *= A68G_REAL_EPS;
270 if (acc < c[1]) {
271 diagnostic (A68G_MATH_WARNING, A68G (f_entry), WARNING_MATH_ACCURACY, NULL);
272 }
273 INT_T N = a68g_round (c[0]);
274 REAL_T err = 0, z = 2 * x, u = 0, v = 0, w = 0;
275 for (int i = 1; i <= N; i++) {
276 if (err > acc) {
277 w = v;
278 v = u;
279 u = z * v - w + c[i];
280 }
281 err += a68g_abs_real (c[i]);
282 }
283 return 0.5 * (u - w);
284 }
285
286 // Compute ln (1 + x) accurately.
287 // Some C99 platforms implement this incorrectly.
288
289 REAL_T a68g_ln1p_real (REAL_T x)
290 {
291 // Based on GNU GSL's gsl_sf_log_1plusx_e.
292 A68G_INVALID (x <= -1);
293 if (a68g_abs_real (x) < pow (A68G_REAL_EPS, 1 / 6.0)) {
294 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;
295 const REAL_T t = c5 + x * (c6 + x * (c7 + x * (c8 + x * c9)));
296 return x * (1 + x * (c1 + x * (c2 + x * (c3 + x * (c4 + x * t)))));
297 } else if (a68g_abs_real (x) < 0.5) {
298 REAL_T t = (8 * x + 1) / (x + 2) / 2;
299 return x * a68g_chebyshev_real (t, c_ln1p, 0.1);
300 } else {
301 return ln (1 + x);
302 }
303 }
304
305 // Compute ln (x), if possible accurately when x ~ 1.
306
307 REAL_T a68g_ln_real (REAL_T x)
308 {
309 A68G_INVALID (x <= 0);
310 #if (A68G_LEVEL >= 3)
311 if (a68g_abs_real (x - 1) < 0.375) {
312 // Double precision x-1 mitigates cancellation error.
313 return a68g_ln1p_real (DBLEQ (x) - 1.0q);
314 } else {
315 return ln (x);
316 }
317 #else
318 return ln (x);
319 #endif
320 }
321
322 // PROC (REAL) REAL exp
323
324 REAL_T a68g_exp_real (REAL_T x)
325 {
326 return exp (x);
327 }
328
329 // OP ** = (REAL, REAL) REAL
330
331 REAL_T a68g_x_up_y (REAL_T x, REAL_T y)
332 {
333 if (y == 0.0) {
334 return 1;
335 } else {
336 return a68g_exp_real (y * a68g_ln_real (x));
337 }
338 }
339
340 // PROC (REAL) REAL csc
341
342 REAL_T a68g_csc_real (REAL_T x)
343 {
344 REAL_T z = sin (x);
345 A68G_OVERFLOW (z == 0);
346 return 1 / z;
347 }
348
349 // PROC (REAL) REAL acsc
350
351 REAL_T a68g_acsc_real (REAL_T x)
352 {
353 A68G_OVERFLOW (x == 0);
354 return asin (1 / x);
355 }
356
357 // PROC (REAL) REAL sec
358
359 REAL_T a68g_sec_real (REAL_T x)
360 {
361 REAL_T z = cos (x);
362 A68G_OVERFLOW (z == 0);
363 return 1 / z;
364 }
365
366 // PROC (REAL) REAL asec
367
368 REAL_T a68g_asec_real (REAL_T x)
369 {
370 A68G_OVERFLOW (x == 0);
371 return acos (1 / x);
372 }
373
374 // PROC (REAL) REAL cot
375
376 REAL_T a68g_cot_real (REAL_T x)
377 {
378 REAL_T z = sin (x);
379 A68G_OVERFLOW (z == 0);
380 return cos (x) / z;
381 }
382
383 // PROC (REAL) REAL acot
384
385 REAL_T a68g_acot_real (REAL_T x)
386 {
387 A68G_OVERFLOW (x == 0);
388 return atan (1 / x);
389 }
390
391 // PROC atan2 (REAL, REAL) REAL
392
393 REAL_T a68g_atan2_real (REAL_T x, REAL_T y)
394 {
395 if (x == 0) {
396 A68G_INVALID (y == 0);
397 return (y > 0 ? CONST_PI_2 : -CONST_PI_2);
398 } else {
399 REAL_T z = atan (ABS (y / x));
400 if (x < 0) {
401 z = CONST_PI - z;
402 }
403 return (y >= 0 ? z : -z);
404 }
405 }
406
407 //! brief PROC (REAL) REAL cas
408
409 REAL_T a68g_cas_real (REAL_T x)
410 {
411 // Hartley kernel, which Hartley named cas (cosine and sine).
412 return cos (x) + sin (x);
413 }
414
415 //! brief PROC (REAL) REAL sindg
416
417 REAL_T a68g_sindg_real (REAL_T x)
418 {
419 return sin (x * CONST_PI_OVER_180);
420 }
421
422 //! brief PROC (REAL) REAL cosdg
423
424 REAL_T a68g_cosdg_real (REAL_T x)
425 {
426 return cos (x * CONST_PI_OVER_180);
427 }
428
429 //! brief PROC (REAL) REAL tandg
430
431 REAL_T a68g_tandg_real (REAL_T x)
432 {
433 return tan (x * CONST_PI_OVER_180);
434 }
435
436 //! brief PROC (REAL) REAL asindg
437
438 REAL_T a68g_asindg_real (REAL_T x)
439 {
440 return asin (x) * CONST_180_OVER_PI;
441 }
442
443 //! brief PROC (REAL) REAL acosdg
444
445 REAL_T a68g_acosdg_real (REAL_T x)
446 {
447 return acos (x) * CONST_180_OVER_PI;
448 }
449
450 //! brief PROC (REAL) REAL atandg
451
452 REAL_T a68g_atandg_real (REAL_T x)
453 {
454 return atan (x) * CONST_180_OVER_PI;
455 }
456
457 // PROC (REAL) REAL cscd
458
459 REAL_T a68g_cscdg_real (REAL_T x)
460 {
461 REAL_T z = a68g_sindg_real (x);
462 A68G_OVERFLOW (z == 0);
463 return 1 / z;
464 }
465
466 // PROC (REAL) REAL acscdg
467
468 REAL_T a68g_acscdg_real (REAL_T x)
469 {
470 A68G_OVERFLOW (x == 0);
471 return a68g_asindg_real (1 / x);
472 }
473
474 // PROC (REAL) REAL secdg
475
476 REAL_T a68g_secdg_real (REAL_T x)
477 {
478 REAL_T z = a68g_cosdg_real (x);
479 A68G_OVERFLOW (z == 0);
480 return 1 / z;
481 }
482
483 // PROC (REAL) REAL asecdg
484
485 REAL_T a68g_asecdg_real (REAL_T x)
486 {
487 A68G_OVERFLOW (x == 0);
488 return a68g_acosdg_real (1 / x);
489 }
490
491 // PROC (REAL) REAL cotdg
492
493 REAL_T a68g_cot_realdg_real (REAL_T x)
494 {
495 REAL_T z = a68g_sindg_real (x);
496 A68G_OVERFLOW (z == 0);
497 return a68g_cosdg_real (x) / z;
498 }
499
500 // PROC (REAL) REAL acotdg
501
502 REAL_T a68g_acotdg_real (REAL_T z)
503 {
504 A68G_OVERFLOW (z == 0);
505 return a68g_atandg_real (1 / z);
506 }
507
508 // @brief PROC (REAL) REAL sinpi
509
510 REAL_T a68g_sinpi_real (REAL_T x)
511 {
512 x = fmod (x, 2);
513 if (x <= -1) {
514 x += 2;
515 } else if (x > 1) {
516 x -= 2;
517 }
518 // x in <-1, 1].
519 if (x == 0 || x == 1) {
520 return 0;
521 } else if (x == 0.5) {
522 return 1;
523 }
524 if (x == -0.5) {
525 return -1;
526 } else {
527 return sin (CONST_PI * x);
528 }
529 }
530
531 // @brief PROC (REAL) REAL cospi
532
533 REAL_T a68g_cospi_real (REAL_T x)
534 {
535 x = fmod (fabs (x), 2);
536 // x in [0, 2>.
537 if (x == 0.5 || x == 1.5) {
538 return 0;
539 } else if (x == 0) {
540 return 1;
541 } else if (x == 1) {
542 return -1;
543 } else {
544 return cos (CONST_PI * x);
545 }
546 }
547
548 // @brief PROC (REAL) REAL tanpi
549
550 REAL_T a68g_tanpi_real (REAL_T x)
551 {
552 x = fmod (x, 1);
553 if (x <= -0.5) {
554 x += 1;
555 } else if (x > 0.5) {
556 x -= 1;
557 }
558 // x in <-1/2, 1/2].
559 A68G_OVERFLOW (x == 0.5);
560 if (x == -0.25) {
561 return -1;
562 } else if (x == 0) {
563 return 0;
564 } else if (x == 0.25) {
565 return 1;
566 } else {
567 return a68g_sinpi_real (x) / a68g_cospi_real (x);
568 }
569 }
570
571 // @brief PROC (REAL) REAL cotpi
572
573 REAL_T a68g_cot_realpi (REAL_T x)
574 {
575 x = fmod (x, 1);
576 if (x <= -0.5) {
577 x += 1;
578 } else if (x > 0.5) {
579 x -= 1;
580 }
581 // x in <-1/2, 1/2].
582 A68G_OVERFLOW (x == 0);
583 if (x == -0.25) {
584 return -1;
585 } else if (x == 0.25) {
586 return 1;
587 } else if (x == 0.5) {
588 return 0;
589 } else {
590 return a68g_cospi_real (x) / a68g_sinpi_real (x);
591 }
592 }
593
594 // @brief PROC (REAL) REAL asinh
595
596 REAL_T a68g_asinh_real (REAL_T x)
597 {
598 REAL_T a = ABS (x), s = (x < 0 ? -1.0 : 1);
599 if (a > 1 / sqrt (A68G_REAL_EPS)) {
600 return (s * (a68g_ln_real (a) + a68g_ln_real (2)));
601 } else if (a > 2) {
602 return (s * a68g_ln_real (2 * a + 1 / (a + sqrt (a * a + 1))));
603 } else if (a > sqrt (A68G_REAL_EPS)) {
604 REAL_T a2 = a * a;
605 return (s * a68g_ln1p_real (a + a2 / (1 + sqrt (1 + a2))));
606 } else {
607 return (x);
608 }
609 }
610
611 // @brief PROC (REAL) REAL acosh
612
613 REAL_T a68g_acosh_real (REAL_T x)
614 {
615 if (x > 1 / sqrt (A68G_REAL_EPS)) {
616 return (a68g_ln_real (x) + a68g_ln_real (2));
617 } else if (x > 2) {
618 return (a68g_ln_real (2 * x - 1 / (sqrt (x * x - 1) + x)));
619 } else if (x > 1) {
620 REAL_T t = x - 1;
621 return (a68g_ln1p_real (t + sqrt (2 * t + t * t)));
622 } else if (x == 1) {
623 return (0);
624 } else {
625 A68G_INVALID (A68G_TRUE);
626 }
627 }
628
629 // @brief PROC (REAL) REAL atanh
630
631 REAL_T a68g_atanh_real (REAL_T x)
632 {
633 REAL_T a = ABS (x);
634 A68G_INVALID (a >= 1);
635 REAL_T s = (x < 0 ? -1 : 1);
636 if (a >= 0.5) {
637 return (s * 0.5 * a68g_ln1p_real (2 * a / (1 - a)));
638 } else if (a > A68G_REAL_EPS) {
639 return (s * 0.5 * a68g_ln1p_real (2 * a + 2 * a * a / (1 - a)));
640 } else {
641 return (x);
642 }
643 }
644
645 //! @brief Inverse complementary error function.
646
647 REAL_T a68g_inverfc_real (REAL_T y)
648 {
649 A68G_INVALID (y < 0 || y > 2);
650 if (y == 0) {
651 return A68G_REAL_MAX;
652 } else if (y == 1) {
653 return 0;
654 } else if (y == 2) {
655 return -A68G_REAL_MAX;
656 } else {
657 // Next is based on code that originally contained following statement:
658 // Copyright (c) 1996 Takuya Ooura. You may use, copy, modify this
659 // code for any purpose and without fee.
660 REAL_T s, t, u, v, x, z;
661 z = (y <= 1 ? y : 2 - y);
662 v = c_inverfc[0] - a68g_ln_real (z);
663 u = sqrt (v);
664 s = (a68g_ln_real (u) + c_inverfc[1]) / v;
665 t = 1 / (u + c_inverfc[2]);
666 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;
667 t = c_inverfc[9] / (x + c_inverfc[9]);
668 u = t - 0.5;
669 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];
670 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]);
671 x += s * (x * s + 1);
672 return (y <= 1 ? x : -x);
673 }
674 }
675
676 //! @brief Inverse error function.
677
678 REAL_T a68g_inverf_real (REAL_T y)
679 {
680 return a68g_inverfc_real (1 - y);
681 }
682
683 //! @brief PROC (REAL, REAL) REAL ln beta
684
685 REAL_T a68g_ln_beta_real (REAL_T a, REAL_T b)
686 {
687 return lgamma (a) + lgamma (b) - lgamma (a + b);
688 }
689
690 //! @brief PROC (REAL, REAL) REAL beta
691
692 REAL_T a68g_beta_real (REAL_T a, REAL_T b)
693 {
694 return a68g_exp_real (a68g_ln_beta_real (a, b));
695 }
696
697 //! brief PROC (INT) REAL fact
698
699 REAL_T a68g_fact_real (INT_T n)
700 {
701 A68G_INVALID (n < 0 || n > A68G_MAX_FAC);
702 return factable[n];
703 }
704
705 //! brief PROC (INT) REAL ln fact
706
707 REAL_T a68g_ln_fact_real (INT_T n)
708 {
709 A68G_INVALID (n < 0);
710 if (n <= A68G_MAX_FAC) {
711 return ln_factable[n];
712 } else {
713 return lgamma (n + 1);
714 }
715 }
716
717 //! @brief PROC choose = (INT n, m) REAL
718
719 REAL_T a68g_choose_real (INT_T n, INT_T m)
720 {
721 A68G_INVALID (n < m);
722 return factable[n] / (factable[m] * factable[n - m]);
723 }
724
725 //! @brief PROC ln choose = (INT n, m) REAL
726
727 REAL_T a68g_ln_choose_real (INT_T n, INT_T m)
728 {
729 A68G_INVALID (n < m);
730 return a68g_ln_fact_real (n) - (a68g_ln_fact_real (m) + a68g_ln_fact_real (n - m));
731 }
732
733 REAL_T a68g_beta_inc_real (REAL_T s, REAL_T t, REAL_T x)
734 {
735 // Incomplete beta function I{x}(s, t).
736 // Continued fraction, see dlmf.nist.gov/8.17; Lentz's algorithm.
737 if (x < 0 || x > 1) {
738 errno = ERANGE;
739 return -1;
740 } else {
741 const INT_T lim = 16 * sizeof (REAL_T);
742 BOOL_T cont = A68G_TRUE;
743 // Rapid convergence when x <= (s+1)/(s+t+2) or else recursion.
744 if (x > (s + 1) / (s + t + 2)) {
745 // B{x}(s, t) = 1 - B{1-x}(t, s)
746 return 1 - a68g_beta_inc_real (s, t, 1 - x);
747 }
748 // Lentz's algorithm for continued fraction.
749 REAL_T W = 1, F = 1, c = 1, d = 0;
750 INT_T N, m;
751 for (N = 0, m = 0; cont && N < lim; N++) {
752 REAL_T T;
753 if (N == 0) {
754 T = 1;
755 } else if (N % 2 == 0) {
756 // d{2m} := x m(t-m)/((s+2m-1)(s+2m))
757 T = x * m * (t - m) / (s + 2 * m - 1) / (s + 2 * m);
758 } else {
759 // d{2m+1} := -x (s+m)(s+t+m)/((s+2m+1)(s+2m))
760 T = -x * (s + m) * (s + t + m) / (s + 2 * m + 1) / (s + 2 * m);
761 m++;
762 }
763 d = 1 / (T * d + 1);
764 c = T / c + 1;
765 F *= c * d;
766 if (F == W) {
767 cont = A68G_FALSE;
768 } else {
769 W = F;
770 }
771 }
772 // I{x}(s,t)=x^s(1-x)^t / s / B(s,t) F
773 REAL_T beta = a68g_exp_real (lgamma (s) + lgamma (t) - lgamma (s + t));
774 return a68g_x_up_y (x, s) * a68g_x_up_y (1 - x, t) / s / beta * (F - 1);
775 }
776 }
|
© 2002-2025 J.M. van der Veer (jmvdveer@xs4all.nl)
|
