## single-math.c

You can download the current version of Algol 68 Genie and its documentation here.

```   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-2023 J. Marcel van der Veer .
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
20 //! License along with this program. If not, see .
21
22 //! @section Synopsis
23 //!
24 //! REAL math stuff supplementing libc.
25
26 // References:
27 //
28 //   Milton Abramowitz and Irene Stegun, Handbook of Mathematical Functions,
29 //   Dover Publications, New York [1970]
30 //   https://en.wikipedia.org/wiki/Abramowitz_and_Stegun
31
32 #include "a68g.h"
33 #include "a68g-genie.h"
34 #include "a68g-prelude.h"
35 #include "a68g-double.h"
36 #include "a68g-numbers.h"
37 #include "a68g-math.h"
38
39 inline REAL_T a68_max (REAL_T x, REAL_T y)
40 {
41   return (x > y ? x : y);
42 }
43
44 inline REAL_T a68_min (REAL_T x, REAL_T y)
45 {
46   return (x < y ? x : y);
47 }
48
49 inline REAL_T a68_sign (REAL_T x)
50 {
51   return (x == 0 ? 0 : (x > 0 ? 1 : -1));
52 }
53
54 inline REAL_T a68_int (REAL_T x)
55 {
56   return (x >= 0 ? (INT_T) x : -(INT_T) - x);
57 }
58
59 inline INT_T a68_round (REAL_T x)
60 {
61   return (INT_T) (x >= 0 ? x + 0.5 : x - 0.5);
62 }
63
64 #define IS_INTEGER(n) (n == a68_int (n))
65
66 inline REAL_T a68_abs (REAL_T x)
67 {
68   return (x >= 0 ? x : -x);
69 }
70
71 REAL_T a68_fdiv (REAL_T x, REAL_T y)
72 {
73 // This is for generating +-INF.
74   return x / y;
75 }
76
77 REAL_T a68_nan (void)
78 {
79   return a68_fdiv (0.0, 0.0);
80 }
81
82 REAL_T a68_posinf (void)
83 {
84   return a68_fdiv (+1.0, 0.0);
85 }
86
87 REAL_T a68_neginf (void)
88 {
89   return a68_fdiv (-1.0, 0.0);
90 }
91
92 // REAL infinity
93
94 void genie_infinity_real (NODE_T * p)
95 {
96   PUSH_VALUE (p, a68_posinf (), A68_REAL);
97 }
98
99 // REAL minus infinity
100
101 void genie_minus_infinity_real (NODE_T * p)
102 {
103   PUSH_VALUE (p, a68_neginf (), A68_REAL);
104 }
105
106 int a68_finite (REAL_T x)
107 {
108 #if defined (HAVE_ISFINITE)
109   return isfinite (x);
110 #elif defined (HAVE_FINITE)
111   return finite (x);
112 #else
113   (void) x;
114   return A68_TRUE;
115 #endif
116 }
117
118 int a68_isnan (REAL_T x)
119 {
120 #if defined (HAVE_ISNAN)
121   return isnan (x);
122 #elif defined (HAVE_IEEE_COMPARISONS)
123   int status = (x != x);
124   return status;
125 #else
126   return A68_FALSE;
127 #endif
128 }
129
130 int a68_isinf (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 (x) && !a68_isnan (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_T x, INT_T n)
218 {
219 // Only positive n.
220   if (n < 0) {
221     return 1 / a68_x_up_n (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 (!finite (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_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_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 *= MATH_EPSILON;
278   if (acc < c[1]) {
279     diagnostic (A68_MATH_WARNING, A68 (f_entry), WARNING_MATH_ACCURACY, NULL);
280   }
281   INT_T i, N = a68_round (c[0]);
282   REAL_T err = 0, z = 2 * x, u = 0, v = 0, w = 0;
283   for (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 (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_T x)
298 {
299 // Based on GNU GSL's gsl_sf_log_1plusx_e.
300   A68_INVALID (x <= -1);
301   if (a68_abs (x) < pow (DBL_EPSILON, 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 (x) < 0.5) {
306     REAL_T t = (8 * x + 1) / (x + 2) / 2;
307     return x * a68_chebyshev (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_T x)
316 {
317   A68_INVALID (x <= 0);
318 #if (A68_LEVEL >= 3)
319   if (a68_abs (x - 1) < 0.375) {
320 // Double precision x-1 mitigates cancellation error.
321     return a68_ln1p (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_T x)
333 {
334   A68_INVALID (x < LOG_DBL_MIN || x > LOG_DBL_MAX);
335   return exp (x);
336 }
337
338 // OP ** = (REAL, REAL) REAL
339
340 REAL_T a68_x_up_y (REAL_T x, REAL_T y)
341 {
342   return a68_exp (y * a68_ln (x));
343 }
344
345 // PROC (REAL) REAL csc
346
347 REAL_T a68_csc (REAL_T x)
348 {
349   REAL_T z = sin (x);
350   A68_OVERFLOW (z == 0);
351   return 1 / z;
352 }
353
354 // PROC (REAL) REAL acsc
355
356 REAL_T a68_acsc (REAL_T x)
357 {
358   A68_OVERFLOW (x == 0);
359   return asin (1 / x);
360 }
361
362 // PROC (REAL) REAL sec
363
364 REAL_T a68_sec (REAL_T x)
365 {
366   REAL_T z = cos (x);
367   A68_OVERFLOW (z == 0);
368   return 1 / z;
369 }
370
371 // PROC (REAL) REAL asec
372
373 REAL_T a68_asec (REAL_T x)
374 {
375   A68_OVERFLOW (x == 0);
376   return acos (1 / x);
377 }
378
379 // PROC (REAL) REAL cot
380
381 REAL_T a68_cot (REAL_T x)
382 {
383   REAL_T z = sin (x);
384   A68_OVERFLOW (z == 0);
385   return cos (x) / z;
386 }
387
388 // PROC (REAL) REAL acot
389
390 REAL_T a68_acot (REAL_T x)
391 {
392   A68_OVERFLOW (x == 0);
393   return atan (1 / x);
394 }
395
396 // PROC atan2 (REAL, REAL) REAL
397
398 REAL_T a68_atan2 (REAL_T x, REAL_T y)
399 {
400   if (x == 0) {
401     A68_INVALID (y == 0);
402     return (y > 0 ? CONST_PI_2 : -CONST_PI_2);
403   } else {
404     REAL_T z = atan (ABS (y / x));
405     if (x < 0) {
406       z = CONST_PI - z;
407     }
408     return (y >= 0 ? z : -z);
409   }
410 }
411
412 //! brief PROC (REAL) REAL sindg
413
414 REAL_T a68_sindg (REAL_T x)
415 {
416   return sin (x * CONST_PI_OVER_180);
417 }
418
419 //! brief PROC (REAL) REAL cosdg
420
421 REAL_T a68_cosdg (REAL_T x)
422 {
423   return cos (x * CONST_PI_OVER_180);
424 }
425
426 //! brief PROC (REAL) REAL tandg
427
428 REAL_T a68_tandg (REAL_T x)
429 {
430   return tan (x * CONST_PI_OVER_180);
431 }
432
433 //! brief PROC (REAL) REAL asindg
434
435 REAL_T a68_asindg (REAL_T x)
436 {
437   return asin (x) * CONST_180_OVER_PI;
438 }
439
440 //! brief PROC (REAL) REAL acosdg
441
442 REAL_T a68_acosdg (REAL_T x)
443 {
444   return acos (x) * CONST_180_OVER_PI;
445 }
446
447 //! brief PROC (REAL) REAL atandg
448
449 REAL_T a68_atandg (REAL_T x)
450 {
451   return atan (x) * CONST_180_OVER_PI;
452 }
453
454 // PROC (REAL) REAL cotdg
455
456 REAL_T a68_cotdg (REAL_T x)
457 {
458   REAL_T z = a68_sindg (x);
459   A68_OVERFLOW (z == 0);
460   return a68_cosdg (x) / z;
461 }
462
463 // PROC (REAL) REAL acotdg
464
465 REAL_T a68_acotdg (REAL_T z)
466 {
467   A68_OVERFLOW (z == 0);
468   return a68_atandg (1 / z);
469 }
470
471 // @brief PROC (REAL) REAL sinpi
472
473 REAL_T a68_sinpi (REAL_T x)
474 {
475   x = fmod (x, 2);
476   if (x <= -1) {
477     x += 2;
478   } else if (x > 1) {
479     x -= 2;
480   }
481 // x in <-1, 1].
482   if (x == 0 || x == 1) {
483     return 0;
484   } else if (x == 0.5) {
485     return 1;
486   }
487   if (x == -0.5) {
488     return -1;
489   } else {
490     return sin (CONST_PI * x);
491   }
492 }
493
494 // @brief PROC (REAL) REAL cospi
495
496 REAL_T a68_cospi (REAL_T x)
497 {
498   x = fmod (fabs (x), 2);
499 // x in [0, 2>.
500   if (x == 0.5 || x == 1.5) {
501     return 0;
502   } else if (x == 0) {
503     return 1;
504   } else if (x == 1) {
505     return -1;
506   } else {
507     return cos (CONST_PI * x);
508   }
509 }
510
511 // @brief PROC (REAL) REAL tanpi
512
513 REAL_T a68_tanpi (REAL_T x)
514 {
515   x = fmod (x, 1);
516   if (x <= -0.5) {
517     x += 1;
518   } else if (x > 0.5) {
519     x -= 1;
520   }
521 // x in <-1/2, 1/2].
522   A68_OVERFLOW (x == 0.5);
523   if (x == -0.25) {
524     return -1;
525   } else if (x == 0) {
526     return 0;
527   } else if (x == 0.25) {
528     return 1;
529   } else {
530     return a68_sinpi (x) / a68_cospi (x);
531   }
532 }
533
534 // @brief PROC (REAL) REAL cotpi
535
536 REAL_T a68_cotpi (REAL_T x)
537 {
538   x = fmod (x, 1);
539   if (x <= -0.5) {
540     x += 1;
541   } else if (x > 0.5) {
542     x -= 1;
543   }
544 // x in <-1/2, 1/2].
545   A68_OVERFLOW (x == 0);
546   if (x == -0.25) {
547     return -1;
548   } else if (x == 0.25) {
549     return 1;
550   } else if (x == 0.5) {
551     return 0;
552   } else {
553     return a68_cospi (x) / a68_sinpi (x);
554   }
555 }
556
557 // @brief PROC (REAL) REAL asinh
558
559 REAL_T a68_asinh (REAL_T x)
560 {
561   REAL_T a = ABS (x), s = (x < 0 ? -1.0 : 1);
562   if (a > 1 / sqrt (DBL_EPSILON)) {
563     return (s * (a68_ln (a) + a68_ln (2)));
564   } else if (a > 2) {
565     return (s * a68_ln (2 * a + 1 / (a + sqrt (a * a + 1))));
566   } else if (a > sqrt (DBL_EPSILON)) {
567     REAL_T a2 = a * a;
568     return (s * a68_ln1p (a + a2 / (1 + sqrt (1 + a2))));
569   } else {
570     return (x);
571   }
572 }
573
574 // @brief PROC (REAL) REAL acosh
575
576 REAL_T a68_acosh (REAL_T x)
577 {
578   if (x > 1 / sqrt (DBL_EPSILON)) {
579     return (a68_ln (x) + a68_ln (2));
580   } else if (x > 2) {
581     return (a68_ln (2 * x - 1 / (sqrt (x * x - 1) + x)));
582   } else if (x > 1) {
583     REAL_T t = x - 1;
584     return (a68_ln1p (t + sqrt (2 * t + t * t)));
585   } else if (x == 1) {
586     return (0);
587   } else {
588     A68_INVALID (A68_TRUE);
589   }
590 }
591
592 // @brief PROC (REAL) REAL atanh
593
594 REAL_T a68_atanh (REAL_T x)
595 {
596   REAL_T a = ABS (x);
597   A68_INVALID (a >= 1);
598   REAL_T s = (x < 0 ? -1 : 1);
599   if (a >= 0.5) {
600     return (s * 0.5 * a68_ln1p (2 * a / (1 - a)));
601   } else if (a > DBL_EPSILON) {
602     return (s * 0.5 * a68_ln1p (2 * a + 2 * a * a / (1 - a)));
603   } else {
604     return (x);
605   }
606 }
607
608 //! @brief Inverse complementary error function.
609
610 REAL_T a68_inverfc (REAL_T y)
611 {
612   A68_INVALID (y < 0 || y > 2);
613   if (y == 0) {
614     return DBL_MAX;
615   } else if (y == 1) {
616     return 0;
617   } else if (y == 2) {
618     return -DBL_MAX;
619   } else {
620 // Next is based on code that originally contained following statement:
621 //   Copyright (c) 1996 Takuya Ooura. You may use, copy, modify this
622 //   code for any purpose and without fee.
623     REAL_T s, t, u, v, x, z;
624     z = (y <= 1 ? y : 2 - y);
625     v = c_inverfc[0] - a68_ln (z);
626     u = sqrt (v);
627     s = (a68_ln (u) + c_inverfc[1]) / v;
628     t = 1 / (u + c_inverfc[2]);
629     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;
630     t = c_inverfc[9] / (x + c_inverfc[9]);
631     u = t - 0.5;
632     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];
633     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 (x * x - c_inverfc[33]);
634     x += s * (x * s + 1);
635     return (y <= 1 ? x : -x);
636   }
637 }
638
639 //! @brief Inverse error function.
640
641 REAL_T a68_inverf (REAL_T y)
642 {
643   return a68_inverfc (1 - y);
644 }
645
646 //! @brief PROC (REAL, REAL) REAL ln beta
647
648 REAL_T a68_ln_beta (REAL_T a, REAL_T b)
649 {
650   return lgamma (a) + lgamma (b) - lgamma (a + b);
651 }
652
653 //! @brief PROC (REAL, REAL) REAL beta
654
655 REAL_T a68_beta (REAL_T a, REAL_T b)
656 {
657   return a68_exp (a68_ln_beta (a, b));
658 }
659
660 //! brief PROC (INT) REAL fact
661
662 REAL_T a68_fact (INT_T n)
663 {
664   A68_INVALID (n < 0 || n > A68_MAX_FAC);
665   return factable[n];
666 }
667
668 //! brief PROC (INT) REAL ln fact
669
670 REAL_T a68_ln_fact (INT_T n)
671 {
672   A68_INVALID (n < 0);
673   if (n <= A68_MAX_FAC) {
674     return ln_factable[n];
675   } else {
676     return lgamma (n + 1);
677   }
678 }
679
680 //! @brief PROC choose = (INT n, m) REAL
681
682 REAL_T a68_choose (INT_T n, INT_T m)
683 {
684   A68_INVALID (n < m);
685   return factable[n] / (factable[m] * factable[n - m]);
686 }
687
688 //! @brief PROC ln choose = (INT n, m) REAL
689
690 REAL_T a68_ln_choose (INT_T n, INT_T m)
691 {
692   A68_INVALID (n < m);
693   return a68_ln_fact (n) - (a68_ln_fact (m) + a68_ln_fact (n - m));
694 }
695
696 REAL_T a68_beta_inc (REAL_T s, REAL_T t, REAL_T x)
697 {
698 // Incomplete beta function I{x}(s, t).
699 // Continued fraction, see dlmf.nist.gov/8.17; Lentz's algorithm.
700   if (x < 0 || x > 1) {
701     errno = ERANGE;
702     return -1;
703   } else {
704     const INT_T lim = 16 * sizeof (REAL_T);
705     BOOL_T cont = A68_TRUE;
706 // Rapid convergence when x <= (s+1)/(s+t+2) or else recursion.
707     if (x > (s + 1) / (s + t + 2)) {
708 // B{x}(s, t) = 1 - B{1-x}(t, s)
709       return 1 - a68_beta_inc (s, t, 1 - x);
710     }
711 // Lentz's algorithm for continued fraction.
712     REAL_T W = 1, F = 1, c = 1, d = 0;
713     INT_T N, m;
714     for (N = 0, m = 0; cont && N < lim; N++) {
715       REAL_T T;
716       if (N == 0) {
717         T = 1;
718       } else if (N % 2 == 0) {
719 // d{2m} := x m(t-m)/((s+2m-1)(s+2m))
720         T = x * m * (t - m) / (s + 2 * m - 1) / (s + 2 * m);
721       } else {
722 // d{2m+1} := -x (s+m)(s+t+m)/((s+2m+1)(s+2m))
723         T = -x * (s + m) * (s + t + m) / (s + 2 * m + 1) / (s + 2 * m);
724         m++;
725       }
726       d = 1 / (T * d + 1);
727       c = T / c + 1;
728       F *= c * d;
729       if (F == W) {
730         cont = A68_FALSE;
731       } else {
732         W = F;
733       }
734     }
735 // I{x}(s,t)=x^s(1-x)^t / s / B(s,t) F
736     REAL_T beta = a68_exp (lgamma (s) + lgamma (t) - lgamma (s + t));
737     return a68_x_up_y (x, s) * a68_x_up_y (1 - x, t) / s / beta * (F - 1);
738   }
739 }
```

 © 2023 J.M. van der Veer jmvdveer@xs4all.nl