From 373efa085dcea5fe6b4539cd875b6bd8645f16fa Mon Sep 17 00:00:00 2001
From: Oliver Kiddle <okiddle@yahoo.co.uk>
Date: Sun, 13 May 2018 10:02:01 +0200
Subject: Nelson H. F. Beebe: 19597 (rebased 42369): return Inf, NaN etc from
 floating point operations instead of errors to allow non-stop IEEE 754
 arithmetic

---
 Src/math.c | 50 ++++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 46 insertions(+), 4 deletions(-)

(limited to 'Src/math.c')

diff --git a/Src/math.c b/Src/math.c
index c3831602b..cdfe80bb4 100644
--- a/Src/math.c
+++ b/Src/math.c
@@ -578,6 +578,37 @@ int outputradix;
 /**/
 int outputunderscore;
 
+#ifndef HAVE_ISINF
+/**/
+int
+isinf(double x)
+{
+    if ((-1.0 < x) && (x < 1.0))       /* x is small, and thus finite */
+       return (0);
+    else if ((x + x) == x)             /* only true if x == Infinity */
+       return (1);
+    else                               /* must be finite (normal or subnormal), or NaN */
+       return (0);
+}
+#endif
+
+#if !defined(HAVE_ISNAN)
+/**/
+static double
+store(double *x)
+{
+    return (*x);
+}
+
+/**/
+int
+isnan(double x)
+{
+    /* (x != x) should be sufficient, but some compilers incorrectly optimize it away */
+    return (store(&x) != store(&x));
+}
+#endif
+
 /**/
 static int
 zzlex(void)
@@ -791,6 +822,21 @@ zzlex(void)
 	    break;
 	/* Fall through! */
 	default:
+	    if (strcmp(ptr-1, "NaN") == 0) {
+		yyval.type = MN_FLOAT;
+		yyval.u.d = 0.0;
+		yyval.u.d /= yyval.u.d;
+		ptr += 2;
+		return NUM;
+	    }
+	    else if (strcmp(ptr-1, "Inf") == 0) {
+		yyval.type = MN_FLOAT;
+		yyval.u.d = 0.0;
+		yyval.u.d = 1.0 / yyval.u.d;
+		ptr += 2;
+		return NUM;
+	    }
+
 	    if (idigit(*--ptr) || *ptr == '.')
 		return lexconstant();
 	    if (*ptr == '#') {
@@ -1068,10 +1114,6 @@ callmathfunc(char *o)
 static int
 notzero(mnumber a)
 {
-    if ((a.type & MN_INTEGER) ? a.u.l == 0 : a.u.d == 0.0) {
-	zerr("division by zero");
-	return 0;
-    }
     return 1;
 }
 
-- 
cgit 1.4.1