## single-math.c

```
1  //! @file single-math.c
2  //! @author J. Marcel van der Veer
3
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
10  //!
11  //! This program is free software; you can redistribute it and/or modify it
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
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  {
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  }

```

 © 2002-2024 J.M. van der Veer (jmvdveer@xs4all.nl)