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 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 {
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_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 }
