From d9a105dabc32bd90487a3c12d6485366b766ffe0 Mon Sep 17 00:00:00 2001
From: Marius Gerbershagen <marius.gerbershagen@gmail.com>
Date: Tue, 5 Mar 2019 20:56:39 +0100
Subject: [PATCH 01/10] numlib.lsp: make signum return proper values for
 infinity/NaN

---
 src/lsp/numlib.lsp | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/src/lsp/numlib.lsp b/src/lsp/numlib.lsp
index 1358f91d7..17c590f87 100644
--- a/src/lsp/numlib.lsp
+++ b/src/lsp/numlib.lsp
@@ -126,7 +126,15 @@ Returns zero for non-complex numbers."
   "Args: (number)
 Returns a number that represents the sign of NUMBER.  Returns NUMBER If it is
 zero.  Otherwise, returns the value of (/ NUMBER (ABS NUMBER))"
-  (if (zerop x) x (/ x (abs x))))
+  (if (complexp x)
+      (cis (atan (imagpart x) (realpart x)))
+      (let ((result (cond ((> x 0) 1)
+                          ((< x 0) -1)
+                          (t ; x is 0 or NaN
+                           x))))
+        (if (floatp x)
+            (float result x)
+            result))))
 
 (defun cis (x)
   "Args: (radians)

From 046c6b9f32e5b7f13a775b44ce00147b761ed5fe Mon Sep 17 00:00:00 2001
From: Marius Gerbershagen <marius.gerbershagen@gmail.com>
Date: Tue, 5 Mar 2019 20:57:52 +0100
Subject: [PATCH 02/10] number_compare.d: fix number comparison for NaN

All comparisons with NaN should return NIL for consistency with the
IEEE standard.
---
 src/c/numbers/number_compare.d | 21 +++++++++++++++------
 src/cmp/sysfun.lsp             |  8 ++++----
 src/h/external.h               |  7 +++++++
 3 files changed, 26 insertions(+), 10 deletions(-)

diff --git a/src/c/numbers/number_compare.d b/src/c/numbers/number_compare.d
index d43e5d11f..d43a9aa59 100644
--- a/src/c/numbers/number_compare.d
+++ b/src/c/numbers/number_compare.d
@@ -159,9 +159,18 @@ monotonic(int s, int t, int narg, ecl_va_list nums)
     FEwrong_type_nth_arg(@[<], 1, c, @[real]);
   }
   /* INV: type check occurs in ecl_number_compare() */
-  for (c = ecl_va_arg(nums); --narg; c = d) {
+  c = ecl_va_arg(nums);
+#ifdef ECL_IEEE_FP
+  if (ecl_float_nan_p(c))
+    return1(ECL_NIL);
+#endif
+  for (; --narg; c = d) {
     d = ecl_va_arg(nums);
-    if (s*ecl_number_compare(d, c) < t)
+    if (
+#ifdef ECL_IEEE_FP
+        ecl_float_nan_p(d) ||
+#endif
+        s*ecl_number_compare(d, c) < t)
       return1(ECL_NIL);
   }
   return1(ECL_T);
@@ -173,7 +182,7 @@ monotonic(int s, int t, int narg, ecl_va_list nums)
     ecl_va_end(nums);                                           \
     return result; }
 
-cl_object @<= MONOTONIC( 1, 0);
-cl_object @>= MONOTONIC(-1, 0);
-cl_object @<  MONOTONIC( 1, 1);
-cl_object @>  MONOTONIC(-1, 1);
+cl_object @<= MONOTONIC( 1, 0)
+cl_object @>= MONOTONIC(-1, 0)
+cl_object @<  MONOTONIC( 1, 1)
+cl_object @>  MONOTONIC(-1, 1)
diff --git a/src/cmp/sysfun.lsp b/src/cmp/sysfun.lsp
index f0a79d473..702010aec 100644
--- a/src/cmp/sysfun.lsp
+++ b/src/cmp/sysfun.lsp
@@ -385,22 +385,22 @@
 (def-inline /= :always (t t) :bool "!ecl_number_equalp(#0,#1)")
 (def-inline /= :always (fixnum-float fixnum-float) :bool "(#0)!=(#1)")
 
-(def-inline < :always (t t) :bool "ecl_number_compare(#0,#1)<0")
+(def-inline < :always (t t) :bool "ecl_lower(#0,#1)")
 (def-inline < :always (fixnum-float fixnum-float) :bool "(#0)<(#1)")
 (def-inline < :always (fixnum-float fixnum-float fixnum-float) :bool
             "@012;((#0)<(#1) && (#1)<(#2))")
 
-(def-inline > :always (t t) :bool "ecl_number_compare(#0,#1)>0")
+(def-inline > :always (t t) :bool "ecl_greater(#0,#1)")
 (def-inline > :always (fixnum-float fixnum-float) :bool "(#0)>(#1)")
 (def-inline > :always (fixnum-float fixnum-float fixnum-float) :bool
             "@012;((#0)>(#1) && (#1)>(#2))")
 
-(def-inline <= :always (t t) :bool "ecl_number_compare(#0,#1)<=0")
+(def-inline <= :always (t t) :bool "ecl_lowereq(#0,#1)")
 (def-inline <= :always (fixnum-float fixnum-float) :bool "(#0)<=(#1)")
 (def-inline <= :always (fixnum-float fixnum-float fixnum-float) :bool
             "@012;((#0)<=(#1) && (#1)<=(#2))")
 
-(def-inline >= :always (t t) :bool "ecl_number_compare(#0,#1)>=0")
+(def-inline >= :always (t t) :bool "ecl_greatereq(#0,#1)")
 (def-inline >= :always (fixnum-float fixnum-float) :bool "(#0)>=(#1)")
 (def-inline >= :always (fixnum-float fixnum-float fixnum-float) :bool
             "@012;((#0)>=(#1) && (#1)>=(#2))")
diff --git a/src/h/external.h b/src/h/external.h
index c671a2e11..fe440988f 100755
--- a/src/h/external.h
+++ b/src/h/external.h
@@ -1191,10 +1191,17 @@ extern ECL_API cl_object cl_min _ECL_ARGS((cl_narg narg, cl_object min, ...));
 
 extern ECL_API int ecl_number_equalp(cl_object x, cl_object y);
 extern ECL_API int ecl_number_compare(cl_object x, cl_object y);
+#ifdef ECL_IEEE_FP
+#define ecl_lowereq(x,y) (!ecl_float_nan_p(x) && !ecl_float_nan_p(y) && ecl_number_compare((x),(y)) <= 0)
+#define ecl_greatereq(x,y) (!ecl_float_nan_p(x) && !ecl_float_nan_p(y) && ecl_number_compare((x),(y)) >= 0)
+#define ecl_lower(x,y) (!ecl_float_nan_p(x) && !ecl_float_nan_p(y) && ecl_number_compare((x),(y)) < 0)
+#define ecl_greater(x,y) (!ecl_float_nan_p(x) && !ecl_float_nan_p(y) && ecl_number_compare((x),(y)) > 0)
+#else
 #define ecl_lowereq(x,y) (ecl_number_compare((x),(y)) <= 0)
 #define ecl_greatereq(x,y) (ecl_number_compare((x),(y)) >= 0)
 #define ecl_lower(x,y) (ecl_number_compare((x),(y)) < 0)
 #define ecl_greater(x,y) (ecl_number_compare((x),(y)) > 0)
+#endif
 
 /* num_log.c */
 

From a597fd5379aea148edbc85ae2ea720720ab9903a Mon Sep 17 00:00:00 2001
From: Marius Gerbershagen <marius.gerbershagen@gmail.com>
Date: Mon, 11 Mar 2019 18:51:42 +0100
Subject: [PATCH 03/10] eql: fix NaN comparison

Old approach doesn't work reliably on x86_64. See e.g.

> (eql (+ ext:double-float-negative-infinity ext:double-float-positive-infinity) (ext:nan))
NIL
---
 src/c/predicate.d | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/c/predicate.d b/src/c/predicate.d
index a14216845..87eb68f3b 100644
--- a/src/c/predicate.d
+++ b/src/c/predicate.d
@@ -258,7 +258,7 @@ cl_eq(cl_object x, cl_object y)
 #define FLOAT_EQL(name, type)                                           \
   static bool name(type a, type b) {                                    \
     if (a == b) return signbit(a) == signbit(b);                        \
-    if (isnan(a) || isnan(b)) return !memcmp(&a, &b, sizeof(type));   \
+    if (isnan(a) || isnan(b)) return isnan(a) && isnan(b);              \
     return 0;                                                           \
   }
 #endif

From 5c7aecc15f3e43b5cad3058aacf9f234f676c090 Mon Sep 17 00:00:00 2001
From: Marius Gerbershagen <marius.gerbershagen@gmail.com>
Date: Sat, 13 Jul 2019 17:43:10 +0200
Subject: [PATCH 04/10] ieee-fp: fix min/max for NaN

As for the C functions fmin/fmax, we ignore NaNs
---
 src/c/numbers/minmax.d | 12 ++++++++++--
 src/cmp/sysfun.lsp     |  6 ++++--
 2 files changed, 14 insertions(+), 4 deletions(-)

diff --git a/src/c/numbers/minmax.d b/src/c/numbers/minmax.d
index dde08dbb8..a691212ac 100644
--- a/src/c/numbers/minmax.d
+++ b/src/c/numbers/minmax.d
@@ -24,7 +24,11 @@
     }
   } else do {
       cl_object numi = ecl_va_arg(nums);
-      if (ecl_number_compare(max, numi) < 0)
+      if (ecl_lower(max, numi)
+#ifdef ECL_IEEE_FP
+          || ecl_float_nan_p(max)
+#endif
+          )
         max = numi;
     } while (--narg);
   @(return max);
@@ -40,7 +44,11 @@
     }
   } else do {
       cl_object numi = ecl_va_arg(nums);
-      if (ecl_number_compare(min, numi) > 0)
+      if (ecl_greater(min, numi)
+#ifdef ECL_IEEE_FP
+          || ecl_float_nan_p(min)
+#endif
+          )
         min = numi;
     } while (--narg);
   @(return min);
diff --git a/src/cmp/sysfun.lsp b/src/cmp/sysfun.lsp
index 702010aec..3f9bf0239 100644
--- a/src/cmp/sysfun.lsp
+++ b/src/cmp/sysfun.lsp
@@ -405,10 +405,12 @@
 (def-inline >= :always (fixnum-float fixnum-float fixnum-float) :bool
             "@012;((#0)>=(#1) && (#1)>=(#2))")
 
-(def-inline max :always (t t) t "@01;(ecl_number_compare(#0,#1)>=0?#0:#1)")
+#+ieee-floating-point (def-inline max :always (t t) t "@01;((ecl_float_nan_p(#1) || ecl_greatereq(#0,#1))?#0:#1)")
+#-ieee-floating-point (def-inline max :always (t t) t "@01;(ecl_greatereq(#0,#1)?#0:#1)")
 (def-inline max :always (fixnum fixnum) :fixnum "@01;(#0)>=(#1)?#0:#1")
 
-(def-inline min :always (t t) t "@01;(ecl_number_compare(#0,#1)<=0?#0:#1)")
+#+ieee-floating-point (def-inline min :always (t t) t "@01;((ecl_float_nan_p(#1) || ecl_lowereq(#0,#1))?#0:#1)")
+#-ieee-floating-point (def-inline min :always (t t) t "@01;(ecl_lowereq(#0,#1)?#0:#1)")
 (def-inline min :always (fixnum fixnum) :fixnum "@01;(#0)<=(#1)?#0:#1")
 
 ;; file num_log.d

From fa28f0877030f972c19f12b7925cde0ce3423110 Mon Sep 17 00:00:00 2001
From: Marius Gerbershagen <marius.gerbershagen@gmail.com>
Date: Sat, 13 Jul 2019 17:45:38 +0200
Subject: [PATCH 05/10] ieee-fp: fix ext:float-{infinity/nan}-p

ext:float-infinity-p returned true for NaNs and ext:float-nan-p was
slower than necessary.
---
 src/c/num_pred.d | 27 +++++++++++++--------------
 1 file changed, 13 insertions(+), 14 deletions(-)

diff --git a/src/c/num_pred.d b/src/c/num_pred.d
index 1ca72747d..edb51e6b2 100644
--- a/src/c/num_pred.d
+++ b/src/c/num_pred.d
@@ -64,17 +64,16 @@ si_float_infinity_p(cl_object x)
 bool
 ecl_float_nan_p(cl_object x)
 {
-  return !ecl_number_equalp(x,x);
-/*   switch (ecl_t_of(x)) { */
-/*   case t_singlefloat: */
-/*     return !isnan(ecl_single_float(x)); */
-/*   case t_doublefloat: */
-/*     return !isnan(ecl_double_float(x)); */
-/*   case t_longfloat: */
-/*     return !isnan(ecl_long_float(x)); */
-/*   default: */
-/*     return 0; */
-/*   } */
+  switch (ecl_t_of(x)) {
+  case t_singlefloat:
+    return isnan(ecl_single_float(x));
+  case t_doublefloat:
+    return isnan(ecl_double_float(x));
+  case t_longfloat:
+    return isnan(ecl_long_float(x));
+  default:
+    return 0;
+  }
 }
 
 bool
@@ -82,11 +81,11 @@ ecl_float_infinity_p(cl_object x)
 {
   switch (ecl_t_of(x)) {
   case t_singlefloat:
-    return !isfinite(ecl_single_float(x));
+    return isinf(ecl_single_float(x));
   case t_doublefloat:
-    return !isfinite(ecl_double_float(x));
+    return isinf(ecl_double_float(x));
   case t_longfloat:
-    return !isfinite(ecl_long_float(x));
+    return isinf(ecl_long_float(x));
   default:
     return 0;
   }

From ac19a3f0a9c15f9a4c2776519813a81cb54edd48 Mon Sep 17 00:00:00 2001
From: Marius Gerbershagen <marius.gerbershagen@gmail.com>
Date: Sat, 13 Jul 2019 21:41:27 +0200
Subject: [PATCH 06/10] ieee-fp: fix round function for NaN

---
 src/c/numbers/round.d | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/c/numbers/round.d b/src/c/numbers/round.d
index fc3bcc308..f9f83fc07 100644
--- a/src/c/numbers/round.d
+++ b/src/c/numbers/round.d
@@ -53,6 +53,8 @@ round_double(double d)
       }
     }
     return q;
+  } else if (isnan(d)) {
+    return d;
   } else {
     return -round_double(-d);
   }
@@ -70,6 +72,8 @@ round_long_double(long double d)
       }
     }
     return q;
+  } else if (isnan(d)) {
+    return d;
   } else {
     return -round_long_double(-d);
   }

From 327da031cd638a52628bc0b45ca9984fd49807fb Mon Sep 17 00:00:00 2001
From: Marius Gerbershagen <marius.gerbershagen@gmail.com>
Date: Sat, 13 Jul 2019 21:41:58 +0200
Subject: [PATCH 07/10] ieee-fp: fix asin/acos/asinh/acosh/atanh for NaN

---
 src/lsp/numlib.lsp | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/lsp/numlib.lsp b/src/lsp/numlib.lsp
index 17c590f87..fa9b087f0 100644
--- a/src/lsp/numlib.lsp
+++ b/src/lsp/numlib.lsp
@@ -151,7 +151,7 @@ RADIANS) and (SIN RADIANS) respectively."
          (setf ,arg (float ,arg)))
        (typecase ,arg
          (single-float
-          (if ,restriction
+          (if (or ,restriction #+ieee-floating-point (ext:float-nan-p ,arg))
               (ffi::c-inline (,arg) (:float) :float
                              ,(format nil "~af(#0)" name)
                              :one-liner t)
@@ -162,7 +162,7 @@ RADIANS) and (SIN RADIANS) respectively."
               #-complex-float
               (progn ,@gencomplex)))
          (double-float
-          (if ,restriction
+          (if (or ,restriction #+ieee-floating-point (ext:float-nan-p ,arg))
               (ffi::c-inline (,arg) (:double) :double
                              ,(format nil "~a(#0)" name)
                              :one-liner t)
@@ -173,7 +173,7 @@ RADIANS) and (SIN RADIANS) respectively."
               #-complex-float
               (progn ,@gencomplex)))
          (long-float
-          (if ,restriction
+          (if (or ,restriction #+ieee-floating-point (ext:float-nan-p ,arg))
               (ffi::c-inline (,arg) (:long-double) :long-double
                              ,(format nil "~al(#0)" name)
                              :one-liner t)

From 57db4c813c09e293ae44813912b140d1df35d08c Mon Sep 17 00:00:00 2001
From: Marius Gerbershagen <marius.gerbershagen@gmail.com>
Date: Wed, 14 Aug 2019 20:41:39 +0200
Subject: [PATCH 08/10] complex-float: use CMPLX for constructing complex
 floats
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Works better for edge cases such as x + I*∞, which otherwise would
lead to a complex with NaN real part.
---
 src/c/number.d   | 12 ++++++------
 src/h/internal.h |  9 +++++++++
 2 files changed, 15 insertions(+), 6 deletions(-)

diff --git a/src/c/number.d b/src/c/number.d
index 7665349fb..2f779b398 100644
--- a/src/c/number.d
+++ b/src/c/number.d
@@ -539,9 +539,9 @@ ecl_make_complex(cl_object r, cl_object i)
   if (!ECL_REAL_TYPE_P(ti)) { ecl_type_error(@'complex', "imaginary part", i, @'real'); }
   switch((tr > ti) ? tr : ti) {
 #ifdef ECL_COMPLEX_FLOAT
-  case t_longfloat:   return ecl_make_clfloat(ecl_to_long_double(r) + I * ecl_to_long_double(i));
-  case t_doublefloat: return ecl_make_cdfloat(ecl_to_double(r) + I * ecl_to_double(i));
-  case t_singlefloat: return ecl_make_csfloat(ecl_to_float(r) + I * ecl_to_float(i));
+  case t_longfloat:   return ecl_make_clfloat(CMPLXL(ecl_to_long_double(r), ecl_to_long_double(i)));
+  case t_doublefloat: return ecl_make_cdfloat(CMPLX(ecl_to_double(r), ecl_to_double(i)));
+  case t_singlefloat: return ecl_make_csfloat(CMPLXF(ecl_to_float(r), ecl_to_float(i)));
 #else
   case t_singlefloat:
     c = ecl_alloc_object(t_complex);
@@ -599,17 +599,17 @@ ecl_make_complex_float(cl_object r, cl_object i)
   case t_singlefloat:
     if (ti != tr) { ecl_type_error(@'si::complex-float',"imag part", i, @'single-float'); }
     result = ecl_alloc_object(t_csfloat);
-    ecl_csfloat(result) = ecl_single_float(r) + ecl_single_float(i) * I;
+    ecl_csfloat(result) = CMPLXF(ecl_single_float(r), ecl_single_float(i));
     break;
   case t_doublefloat:
     if (ti != tr) { ecl_type_error(@'si::complex-float',"imag part", i, @'double-float'); }
     result = ecl_alloc_object(t_cdfloat);
-    ecl_cdfloat(result) = ecl_double_float(r) + ecl_double_float(i) * I;
+    ecl_cdfloat(result) = CMPLX(ecl_double_float(r), ecl_double_float(i));
     break;
   case t_longfloat:
     if (ti != tr) { ecl_type_error(@'si::complex-float',"imag part", i, @'long-float'); }
     result = ecl_alloc_object(t_clfloat);
-    ecl_clfloat(result) = ecl_long_float(r) + ecl_long_float(i) * I;
+    ecl_clfloat(result) = CMPLXL(ecl_long_float(r), ecl_long_float(i));
     break;
   default:
     ecl_type_error(@'si::complex-float',"real part", r, @'float');
diff --git a/src/h/internal.h b/src/h/internal.h
index f163c6b3f..8892e1190 100755
--- a/src/h/internal.h
+++ b/src/h/internal.h
@@ -679,6 +679,15 @@ static union {
 
 #ifdef ECL_COMPLEX_FLOAT
 #include <complex.h>
+#ifndef CMPLXF
+# define CMPLXF(x, y) ((float complex)((float)(x) + I * (float)(y)))
+#endif
+#ifndef CMPLX
+# define CMPLX(x, y) ((double complex)((double)(x) + I * (double)(y)))
+#endif
+#ifndef CMPLXL
+# define CMPLXL(x, y) ((long double complex)((long double)(x) + I * (long double)(y)))
+#endif
 #endif
 
 #ifdef __cplusplus

From 3c2105fe3d089a59dabec850969308c613a74d00 Mon Sep 17 00:00:00 2001
From: Marius Gerbershagen <marius.gerbershagen@gmail.com>
Date: Wed, 27 Mar 2019 21:56:57 +0100
Subject: [PATCH 09/10] ieee-fp: add tests for numeric functions with
 infinity/NaN

---
 src/c/numbers/ceiling.d            |   2 +-
 src/c/numbers/floor.d              |   2 +-
 src/c/numbers/round.d              |   2 +-
 src/c/numbers/truncate.d           |   2 +-
 src/lsp/numlib.lsp                 |   2 +-
 src/tests/ecl-tests.lisp           |   3 +-
 src/tests/normal-tests/ieee-fp.lsp | 650 ++++++++++++++++++++++++++++-
 src/tests/universe.lisp            |   6 +
 8 files changed, 641 insertions(+), 28 deletions(-)

diff --git a/src/c/numbers/ceiling.d b/src/c/numbers/ceiling.d
index 750b70f6e..4599cca33 100644
--- a/src/c/numbers/ceiling.d
+++ b/src/c/numbers/ceiling.d
@@ -21,7 +21,7 @@
 @(defun ceiling (x &optional (y OBJNULL))
   @
   if (narg == 1)
-  return ecl_ceiling1(x);
+    return ecl_ceiling1(x);
   else
     return ecl_ceiling2(x, y);
   @)
diff --git a/src/c/numbers/floor.d b/src/c/numbers/floor.d
index b2cb828d8..7648473cf 100644
--- a/src/c/numbers/floor.d
+++ b/src/c/numbers/floor.d
@@ -24,7 +24,7 @@
 @(defun floor (x &optional (y OBJNULL))
   @
   if (narg == 1)
-  return ecl_floor1(x);
+    return ecl_floor1(x);
   else
     return ecl_floor2(x, y);
   @)
diff --git a/src/c/numbers/round.d b/src/c/numbers/round.d
index f9f83fc07..d46eba28d 100644
--- a/src/c/numbers/round.d
+++ b/src/c/numbers/round.d
@@ -26,7 +26,7 @@
 @(defun round (x &optional (y OBJNULL))
   @
   if (narg == 1)
-  return ecl_round1(x);
+    return ecl_round1(x);
   else
     return ecl_round2(x, y);
   @)
diff --git a/src/c/numbers/truncate.d b/src/c/numbers/truncate.d
index 65a573c9c..d3a705db9 100644
--- a/src/c/numbers/truncate.d
+++ b/src/c/numbers/truncate.d
@@ -79,7 +79,7 @@ ecl_truncate2(cl_object x, cl_object y)
 @(defun truncate (x &optional (y OBJNULL))
   @
   if (narg == 1)
-  return ecl_truncate1(x);
+    return ecl_truncate1(x);
   else
     return ecl_truncate2(x, y);
   @)
diff --git a/src/lsp/numlib.lsp b/src/lsp/numlib.lsp
index fa9b087f0..32ba1d1fd 100644
--- a/src/lsp/numlib.lsp
+++ b/src/lsp/numlib.lsp
@@ -319,7 +319,7 @@ Returns the hyperbolic arc tangent of NUMBER."
   (declare (number z) (si::c-local))
   (/ (- (log (1+ z)) (log (- 1 z))) 2))
 
-(defun ffloor (x &optional (y 1))
+(defun ffloor (x &optional (y 1.0f0))
   "Args: (number &optional (divisor 1))
 Same as FLOOR, but returns a float as the first value."
   (multiple-value-bind (i r) (floor x y)
diff --git a/src/tests/ecl-tests.lisp b/src/tests/ecl-tests.lisp
index e18755b2e..cccf9751c 100644
--- a/src/tests/ecl-tests.lisp
+++ b/src/tests/ecl-tests.lisp
@@ -121,7 +121,8 @@ as a second value."
 
 ;;; Approximate equality function
 (defun approx= (x y &optional (eps (epsilon x)))
-  (<= (abs (/ (- x y) (max (abs x) 1))) eps))
+  (or (= x y)
+      (<= (abs (/ (- x y) (max (abs x) 1))) eps)))
 
 (defun epsilon (number)
   (etypecase number
diff --git a/src/tests/normal-tests/ieee-fp.lsp b/src/tests/normal-tests/ieee-fp.lsp
index 1b5d58793..401a6187d 100644
--- a/src/tests/normal-tests/ieee-fp.lsp
+++ b/src/tests/normal-tests/ieee-fp.lsp
@@ -18,6 +18,113 @@
             ,@body)
        (si:trap-fpe bits t))))
 
+(defconstant +number-subtypes-and-values+
+  '((fixnum 3)
+    (bignum #.(- most-negative-fixnum 13))
+    (ratio 3/5)
+    (single-float 1.23f0)
+    (double-float 4.56d0)
+    (long-float 7.89l0)
+    ((complex fixnum) #C(2 3))
+    ((complex integer) #.(complex 2 (+ most-positive-fixnum 3)))
+    ((complex rational) #C(5/7 -4))
+    ((complex single-float) #C(0.98f0 77.7f0))
+    ((complex double-float) #C(0.87d0 88.8d0))
+    ((complex long-float) #C(0.76l0 99.9l0))))
+
+(defmacro for-all-number-subtypes ((variable &optional (type number) default-value) &body body)
+  (when (not (subtypep type 'number))
+    (error "~a must be a subtype of number" type))
+  `(progn ,@(loop :for (subtype value) :in +number-subtypes-and-values+
+               :if (subtypep subtype type)
+               :collect `(let ((,variable
+                                ,(if default-value
+                                     `(coerce ,default-value ',subtype)
+                                     value)))
+                           ,@body))))
+
+(defmacro for-all-NaNs (variable &body body)
+  `(loop :for ,variable :in *ieee-fp.nan* :do
+      ,@body))
+
+(defmacro for-all-infinities (minus-infinity-var plus-infinity-var &body body)
+  `(loop :for ,minus-infinity-var :in (remove-if-not #'minusp *ieee-fp.inf*) :do
+      (loop :for ,plus-infinity-var :in (remove-if-not #'plusp *ieee-fp.inf*) :do
+         ,@body)))
+
+(defun type-contagion (&rest numbers)
+  (labels ((type-of-proper (number)
+             (etypecase number
+               (single-float 'single-float)
+               (double-float 'double-float)
+               (long-float 'long-float)
+               (ratio 'ratio)
+               (integer 'integer)
+               (complex `(complex ,(type-of-proper (realpart number))))))
+           (type-contagion1 (t1 &optional t2)
+             (cond ((not t2) t1)
+                   ((subtypep t1 'float)
+                    (cond ((subtypep t2 'float)
+                           (cond ((or (eql t1 'long-float)
+                                      (eql t2 'long-float))
+                                  'long-float)
+                                 ((or (eql t1 'double-float)
+                                      (eql t2 'double-float))
+                                  'double-float)
+                                 ((or (eql t1 'single-float)
+                                      (eql t2 'single-float))
+                                  'single-float)))
+                          ((subtypep t2 'rational)
+                           t1)
+                          ((subtypep t2 'complex)
+                           `(complex ,(type-contagion1 t1 (second t2))))))
+                   ((subtypep t1 'ratio)
+                    (cond ((subtypep t2 'float) t2)
+                          ((subtypep t2 'rational) 'ratio)
+                          ((subtypep t2 'complex)
+                           `(complex ,(type-contagion1 t1 (second t2))))))
+                   ((subtypep t1 'integer)
+                    (cond ((or (subtypep t2 'float) (subtypep t2 'ratio)) t2)
+                          ((subtypep t2 'integer) 'integer)
+                          ((subtypep t2 'complex)
+                           `(complex ,(type-contagion1 t1 (second t2))))))
+                   ((subtypep t1 'complex)
+                    (cond ((or (subtypep t2 'float) (subtypep t2 'rational))
+                           `(complex ,(type-contagion (second t1) t2)))
+                          ((subtypep t2 'complex)
+                           `(complex ,(type-contagion1 (second t1) (second t2)))))))))
+    (reduce #'type-contagion1 (mapcar #'type-of-proper numbers))))
+
+(defmacro type-and-value-check ((function &rest args) &rest values)
+  "Arguments ((function &rest args) &rest values)
+Check that function applied to args returns numbers that are
+mathematically equal to values and have the correct type according to
+Common Lisp type contagion rules."
+  (if (eql function 'funcall)
+      (setf function (first args)
+            args (rest args))
+      (setf function `(quote ,function)))
+  `(is (every #'eql
+              (multiple-value-list (funcall ,function ,@args))
+              (mapcar #'(lambda (v)
+                          (coerce v (type-contagion ,@args)))
+                      (list ,@values)))))
+
+;;; Constructing complex floats with infinities is difficult, because
+;;; any finite real number x plus infinite imaginary number leads to
+;;; an complex number with NaN real part due to a very literal
+;;; application of calculation rules:
+;;; x + I*∞ = x + (0 + I*1)*∞ = x + NaN + I*∞ = NaN + I*∞ 
+;;; Therefore, when comparing the output of numeric functions
+;;; returning complex infinities to expected results, we use the
+;;; following function
+(defmacro complex-equality (equality-comparison form realpart imagpart)
+  `(progn
+     (is (funcall ,equality-comparison (realpart ,form) ,realpart))
+     (is (funcall ,equality-comparison (imagpart ,form) ,imagpart))))
+
+
+
 (test ieee-fp.0001.infinity-eql
   (without-fpe-traps
     (let ((sfni ext:single-float-negative-infinity)
@@ -55,17 +162,280 @@
   (is (ext:float-infinity-p ext:single-float-positive-infinity))
   (is (ext:float-nan-p (si:nan))))
 
-;;; Reported by: Robert Dodier
-;;; URL: https://gitlab.com/embeddable-common-lisp/ecl/issues/299
-(test ieee-fp.0003.b299
-  (finishes (< ext:double-float-negative-infinity 1/3))
-  (finishes (> ext:double-float-negative-infinity 1/3))
-  (finishes (< 1/3 ext:double-float-negative-infinity))
-  (finishes (> 1/3 ext:double-float-negative-infinity))
-  (finishes (< ext:double-float-negative-infinity (1+ most-positive-fixnum)))
-  (finishes (> ext:double-float-negative-infinity (1+ most-positive-fixnum)))
-  (finishes (< (1+ most-positive-fixnum) ext:double-float-negative-infinity))
-  (finishes (> (1+ most-positive-fixnum) ext:double-float-negative-infinity)))
+(test ieee-fp.0004.nan-comparison
+      (for-all-number-subtypes (x real)
+        (for-all-NaNs NaN
+          (is (not (< NaN x)))
+          (is (not (< x NaN)))
+          (is (not (<= NaN x)))
+          (is (not (<= x NaN)))
+          (is (not (> NaN x)))
+          (is (not (> x NaN)))
+          (is (not (>= NaN x)))
+          (is (not (>= x NaN))))))
+
+(test ieee-fp.0005.infinity-comparison
+      (for-all-number-subtypes (x real)
+        (for-all-infinities -infinity +infinity
+          (progn (is (< -infinity x))
+                 (is (not (< x -infinity)))
+                 (is (< x +infinity))
+                 (is (not (< +infinity x))))
+          (progn (is (<= -infinity x))
+                 (is (not (<= x -infinity)))
+                 (is (<= x +infinity))
+                 (is (not (<= +infinity x))))
+          (progn (is (not (> -infinity x)))
+                 (is (> x -infinity))
+                 (is (not (> x +infinity)))
+                 (is (> +infinity x)))
+          (progn (is (not (>= -infinity x)))
+                 (is (>= x -infinity))
+                 (is (not (>= x +infinity)))
+                 (is (>= +infinity x))))))
+
+
+
+;;; We treat NaNs as missing data (like fmin/fmax C functions)
+(test ieee-fp.0006.NaN-min/max
+      (for-all-NaNs NaN
+        (is (eql (min NaN) NaN))
+        (is (eql (max NaN) NaN))
+        (is (eql (min NaN NaN) NaN))
+        (is (eql (max NaN NaN) NaN))
+        (for-all-number-subtypes (x real)
+          (let* ((y (* 2 x))
+                 (min (min x y))
+                 (max (max x y)))
+            (is (= (min x NaN) x))
+            (is (= (max x NaN) x))
+            (is (= (min x y NaN) min))
+            (is (= (max x y NaN) max))))))
+
+(test ieee-fp.0007.infinity-min/max
+      (for-all-number-subtypes (x real)
+        (let* ((y (* 2 x))
+               (max (max x y))
+               (min (min x y)))
+          (for-all-infinities -infinity +infinity
+            (is (= (min -infinity x y) -infinity))
+            (is (= (min +infinity x y) min))
+            (is (= (max -infinity x y) max))
+            (is (= (max +infinity x y) +infinity))))))
+
+(test ieee-fp.0008.NaN-minusp/zerop/plusp
+      (without-fpe-traps
+        (for-all-NaNs NaN
+          (is (not (minusp NaN)))
+          (is (not (zerop NaN)))
+          (is (not (plusp NaN))))))
+
+(test ieee-fp.0009.infinity-minusp/zerop/plusp
+      (for-all-infinities -infinity +infinity
+        (is (minusp -infinity))
+        (is (not (minusp +infinity)))
+        (is (not (zerop -infinity)))
+        (is (not (zerop +infinity)))
+        (is (not (plusp -infinity)))
+        (is (plusp +infinity))))
+
+
+
+;;; If we supported NaN/infinity for rounding functions operating
+;;; on floats, sensible behaviour would be as follows ...
+
+#|
+(test ieee-fp.0010.NaN-ffloor/fceiling/ftruncate/fround
+      (without-fpe-traps
+          (loop :for function :in '(ffloor fceiling ftruncate fround)
+             :do
+             (for-all-number-subtypes (x real)
+               (for-all-NaNs NaN
+                 (type-and-value-check (funcall function NaN) NaN NaN)
+                 (type-and-value-check (funcall function x NaN) NaN NaN)
+                 (type-and-value-check (funcall function NaN x) NaN NaN)
+                 (type-and-value-check (funcall function NaN NaN) NaN NaN)))
+             (for-all-number-subtypes (x float 0)
+               (type-and-value-check (funcall function x x) (ext:nan) (ext:nan))))))
+
+(test ieee-fp.0011.infinity-ffloor/fceiling/ftruncate/fround
+      (without-fpe-traps
+        (loop :for function :in '(ffloor fceiling ftruncate fround)
+           :do
+           (for-all-infinities -infinity +infinity
+             ;; +-infinity and 1
+             (type-and-value-check (funcall function -infinity)
+                                   -infinity (ext:nan))
+             (type-and-value-check (funcall function +infinity)
+                                   +infinity (ext:nan))
+             (for-all-number-subtypes (x real)
+               ;; +-infinity and finite real
+               (type-and-value-check (funcall function -infinity (abs x))
+                                     -infinity (ext:nan))
+               (type-and-value-check (funcall function -infinity (- (abs x)))
+                                     +infinity (ext:nan))
+               (type-and-value-check (funcall function +infinity (abs x))
+                                     +infinity (ext:nan))
+               (type-and-value-check (funcall function +infinity (- (abs x)))
+                                     -infinity (ext:nan))
+               ;; finite real and +-infinity
+               (type-and-value-check (funcall function (abs x) -infinity)
+                                     (case function
+                                       (ffloor -1)
+                                       (t 0))
+                                     (case function
+                                       (ffloor -infinity)
+                                       (t x)))
+               (type-and-value-check (funcall function (- (abs x)) -infinity)
+                                     (case function
+                                       (fceiling 1)
+                                       (t 0))
+                                     (case function
+                                       (fceiling +infinity)
+                                       (t x)))
+               (type-and-value-check (funcall function (abs x) +infinity)
+                                     (case function
+                                       (fceiling 1)
+                                       (t 0))
+                                     (case function
+                                       (fceiling -infinity)
+                                       (t x)))
+               (type-and-value-check (funcall function (- (abs x)) +infinity)
+                                     (case function
+                                       (ffloor -1)
+                                       (t 0))
+                                     (case function
+                                       (ffloor +infinity)
+                                       (t x))))
+             ;; +-infinity and +-infinity
+             (type-and-value-check (funcall function -infinity -infinity)
+                                   (ext:nan) (ext:nan))
+             (type-and-value-check (funcall function -infinity +infinity)
+                                   (ext:nan) (ext:nan))
+             (type-and-value-check (funcall function +infinity -infinity)
+                                   (ext:nan) (ext:nan))
+             (type-and-value-check (funcall function +infinity +infinity)
+                                   (ext:nan) (ext:nan))
+             (for-all-number-subtypes (x float 0)
+               ;; finite real and 0
+               (for-all-number-subtypes (y real)
+                 (type-and-value-check (funcall function (abs y) x)
+                                       +infinity (ext:nan))
+                 (type-and-value-check (funcall function (- (abs y)) x)
+                                       -infinity (ext:nan)))
+               ;; +-infinity and 0
+               (type-and-value-check (funcall function -infinity 0)
+                                     -infinity (ext:nan))
+               (type-and-value-check (funcall function +infinity 0)
+                                     +infinity (ext:nan)))))))
+|#
+
+;;; ... but we don't, therefore everything throws arithmetic errors.
+
+(test ieee-fp.0010.NaN-floor/ceiling/truncate/round/mod/rem
+      (loop :for function :in '(floor ceiling truncate round
+                                ffloor fceiling ftruncate fround
+                                mod rem)
+         :do
+         (for-all-NaNs NaN
+           (unless (member function '(mod rem))
+             (signals arithmetic-error (funcall function NaN)))
+           (signals arithmetic-error (funcall function NaN NaN))
+           (for-all-number-subtypes (x real)
+             (signals arithmetic-error (funcall function x NaN))
+             (signals arithmetic-error (funcall function NaN x))))
+         (for-all-number-subtypes (x float 0)
+           (signals arithmetic-error (funcall function x x)))))
+
+(test ieee-fp.0011.infinity-floor/ceiling/truncate/round
+      (loop :for function :in '(floor ceiling truncate round
+                                ffloor fceiling ftruncate fround
+                                mod rem)
+         :do
+         (for-all-infinities -infinity +infinity
+           ;; +-infinity and 1
+           (unless (member function '(mod rem))
+             (signals arithmetic-error (funcall function -infinity))
+             (signals arithmetic-error (funcall function +infinity)))
+           (for-all-number-subtypes (x real)
+             ;; +-infinity and finite real
+             (signals arithmetic-error (funcall function -infinity x))
+             (signals arithmetic-error (funcall function +infinity x))
+             ;; finite real and +-infinity
+             (signals arithmetic-error (funcall function x -infinity))
+             (signals arithmetic-error (funcall function x +infinity)))
+           ;; +-infinity and +-infinity
+           (signals arithmetic-error (funcall function -infinity -infinity))
+           (signals arithmetic-error (funcall function -infinity +infinity))
+           (signals arithmetic-error (funcall function +infinity -infinity))
+           (signals arithmetic-error (funcall function +infinity +infinity))
+           (for-all-number-subtypes (x float 0)
+             ;; finite real and 0
+             (for-all-number-subtypes (y real)
+               (signals arithmetic-error (funcall function y x)))
+             ;; +-infinity and 0
+             (signals arithmetic-error (funcall function -infinity x))
+             (signals arithmetic-error (funcall function +infinity x))))))
+
+
+
+;;; All numeric functions should return a NaN if given one
+(test ieee-fp.0012.sticky-NaN
+      (without-fpe-traps
+        (for-all-NaNs NaN
+          (macrolet ((f-NaN (&rest functions)
+                       `(progn ,@(loop :for f :in functions :collect
+                                      `(is (eql (,f NaN) NaN))))))
+            (f-NaN sin cos tan
+                   asin acos atan
+                   sinh cosh tanh
+                   asinh acosh atanh
+                   * + - / 1+ 1-
+                   abs signum exp log sqrt
+                   conjugate phase))
+          (for-all-number-subtypes (x real)
+            (macrolet ((f-NaN (&rest functions)
+                         `(progn ,@(loop :for f :in functions :append
+                                        `((type-and-value-check (,f x NaN) NaN)
+                                          (type-and-value-check (,f NaN x) NaN)
+                                          (type-and-value-check (,f NaN NaN) NaN))))))
+              (f-NaN atan * + - /))
+            (macrolet ((f-NaN (&rest functions)
+                         `(progn ,@(loop :for f :in functions :append
+                                        `((type-and-value-check (,f (abs x) NaN) NaN)
+                                          (type-and-value-check (,f NaN (abs x)) NaN)
+                                          (type-and-value-check (,f NaN NaN) NaN))))))
+              (f-NaN expt log)))
+          (is (eql (cis NaN)
+                   (complex (coerce (ext:nan) (type-of NaN))
+                            (coerce (ext:nan) (type-of NaN))))))))
+
+
+
+(test ieee-fp.0013.infinity-sin/cos/tan
+      (without-fpe-traps
+        (for-all-infinities -infinity +infinity
+          (type-and-value-check (sin -infinity) (ext:nan))
+          (type-and-value-check (sin +infinity) (ext:nan))
+          #|(type-and-value-check (sin (complex 0 -infinity)) (complex 0 -infinity))
+          (type-and-value-check (sin (complex 0 +infinity)) (complex 0 +infinity))|#
+          (type-and-value-check (cos -infinity) (ext:nan))
+          (type-and-value-check (cos +infinity) (ext:nan))
+          #|(type-and-value-check (cos (complex 0 -infinity)) +infinity)
+          (type-and-value-check (cos (complex 0 +infinity)) +infinity)|#
+          (type-and-value-check (tan -infinity) (ext:nan))
+          (type-and-value-check (tan +infinity) (ext:nan))
+          #|(type-and-value-check (tan (complex 0 -infinity)) (complex 0 -1))
+          (type-and-value-check (tan (complex 0 +infinity)) (complex 0 +1))|#)))
+
+#|(test ieee-fp.0014.infinity-asin/acos/atan
+      (for-all-infinities -infinity +infinity
+        (complex-equality #'approx= (asin -infinity) (- (/ pi 2)) +infinity)
+        (complex-equality #'approx= (asin +infinity) (/ pi 2) +infinity)
+        (complex-equality #'approx= (acos -infinity) pi -infinity)
+        (complex-equality #'approx= (acos +infinity) 0 -infinity)
+        (is (approx= (atan -infinity) (- (/ pi 2))))
+        (is (approx= (atan +infinity) (/ pi 2)))))|#
 
 
 ;;; Reported by: Raymond Toy
@@ -79,16 +449,7 @@
       (+zerop (elt) (and (zerop elt) (plusp (float-sign elt))))
       (-zerop (elt) (and (zerop elt) (minusp (float-sign elt)))))
 
-  (test ieee-fp.0004.atan2-special-cases.nan-arg
-        (without-fpe-traps
-            ;; (atan (anything) nan) -> nan
-            ;; (atan nan (anything)) -> nan
-            (map () (lambda (n)
-                      (is (si:float-nan-p (atan n (si:nan))))
-                      (is (si:float-nan-p (atan (si:nan) n))))
-                 (append *floats* *ieee-fp.inf* (list (si:nan))))))
-
-  (test ieee-fp.0005.atan2-special-case.zero-arg
+  (test ieee-fp.0015.atan2-special-case.zero-arg
         (+zerop (atan +0.0 +0.0))
         (-zerop (atan -0.0 +0.0))
         (+pi-p  (atan +0.0 -0.0))
@@ -108,7 +469,7 @@
                     (is (-pi/2-p (atan (- n) -0.0))))
                (remove-if-not #'plusp (append *floats* *ieee-fp.inf*)))))
 
-  (test ieee-fp.0006.atan2-special-case.inf-arg
+  (test ieee-fp.0016.atan2-special-case.inf-arg
         ;; (atan +-inf +inf) -> +-pi/4
         (let ((+inf ext:single-float-positive-infinity)
               (-inf ext:single-float-negative-infinity))
@@ -125,3 +486,248 @@
                       (+pi/2-p (atan +inf n))
                       (-pi/2-p (atan -inf n))))
                (remove-if-not #'plusp *floats*)))))
+
+(test ieee-fp.0017.infinity-sinh/cosh/tanh
+      (for-all-infinities -infinity +infinity
+        (type-and-value-check (sinh -infinity) -infinity)
+        (type-and-value-check (sinh +infinity) +infinity)
+        (type-and-value-check (cosh -infinity) +infinity)
+        (type-and-value-check (cosh +infinity) +infinity)
+        (type-and-value-check (tanh -infinity) -1)
+        (type-and-value-check (tanh +infinity) +1)
+        #|(without-fpe-traps
+          (type-and-value-check (sinh (complex 0 -infinity)) (complex (ext:nan) (ext:nan)))
+          (type-and-value-check (sinh (complex 0 +infinity)) (complex (ext:nan) (ext:nan)))
+          (type-and-value-check (cosh (complex 0 -infinity)) (complex (ext:nan) (ext:nan)))
+          (type-and-value-check (cosh (complex 0 +infinity)) (complex (ext:nan) (ext:nan)))
+          (type-and-value-check (tanh (complex 0 -infinity)) (complex (ext:nan) (ext:nan)))
+          (type-and-value-check (tanh (complex 0 +infinity)) (complex (ext:nan) (ext:nan))))|#))
+
+(test ieee-fp.0018.infinity-asinh/acosh/atanh
+      (for-all-infinities -infinity +infinity
+        (type-and-value-check (asinh -infinity) -infinity)
+        (type-and-value-check (asinh +infinity) +infinity)
+        (without-fpe-traps
+          (complex-equality #'approx= (acosh -infinity) +infinity pi))
+        (type-and-value-check (acosh +infinity) +infinity)
+        #| (complex-equality #'approx= (atanh -infinity) 0 (/ pi 2))
+        (complex-equality #'approx= (atanh +infinity) 0 (/ pi 2))|#))
+
+(test ieee-fp.0019.infinity-+/-
+      (for-all-infinities -infinity +infinity
+        (type-and-value-check (+ -infinity) -infinity)
+        (type-and-value-check (+ +infinity) +infinity)
+
+        (type-and-value-check (+ -infinity -infinity) -infinity)
+        (without-fpe-traps
+          (type-and-value-check (+ -infinity +infinity) (ext:nan))
+          (type-and-value-check (+ +infinity -infinity) (ext:nan)))
+        (type-and-value-check (+ +infinity +infinity) +infinity)
+
+        (type-and-value-check (- -infinity) +infinity)
+        (type-and-value-check (- +infinity) -infinity)
+
+        (without-fpe-traps
+          (type-and-value-check (- -infinity -infinity) (ext:nan)))
+        (type-and-value-check (- +infinity -infinity) +infinity)
+        (type-and-value-check (- -infinity +infinity) -infinity)
+        (without-fpe-traps
+          (type-and-value-check (- +infinity +infinity) (ext:nan)))
+        (for-all-number-subtypes (x real)
+          (type-and-value-check (+ x -infinity) -infinity)
+          (type-and-value-check (+ x +infinity) +infinity)
+          (type-and-value-check (+ -infinity x) -infinity)
+          (type-and-value-check (+ +infinity x) +infinity)
+
+          (type-and-value-check (- x -infinity) +infinity)
+          (type-and-value-check (- x +infinity) -infinity)
+          (type-and-value-check (- -infinity x) -infinity)
+          (type-and-value-check (- +infinity x) +infinity))))
+
+(test ieee-fp.0020.infinity-*//
+      (for-all-infinities -infinity +infinity
+        (type-and-value-check (* -infinity) -infinity)
+        (type-and-value-check (* +infinity) +infinity)
+
+        (type-and-value-check (* -infinity -infinity) +infinity)
+        (type-and-value-check (* -infinity +infinity) -infinity)
+        (type-and-value-check (* +infinity -infinity) -infinity)
+        (type-and-value-check (* +infinity +infinity) +infinity)
+
+        (is (eql (/ -infinity) (coerce -0.0 (type-of -infinity))))
+        (is (eql (/ +infinity) (coerce +0.0 (type-of +infinity))))
+
+        (without-fpe-traps
+          (type-and-value-check (/ -infinity -infinity) (ext:nan))
+          (type-and-value-check (/ -infinity +infinity) (ext:nan))
+          (type-and-value-check (/ +infinity -infinity) (ext:nan))
+          (type-and-value-check (/ +infinity +infinity) (ext:nan)))
+        (for-all-number-subtypes (x real)
+          (type-and-value-check (* (abs x) -infinity) -infinity)
+          (type-and-value-check (* (- (abs x)) -infinity) +infinity)
+          (type-and-value-check (* (abs x) +infinity) +infinity)
+          (type-and-value-check (* (- (abs x)) +infinity) -infinity)
+
+          (type-and-value-check (* -infinity (abs x)) -infinity)
+          (type-and-value-check (* -infinity (- (abs x))) +infinity)
+          (type-and-value-check (* +infinity (abs x)) +infinity)
+          (type-and-value-check (* +infinity (- (abs x))) -infinity)
+
+          (is (eql (/ (abs x) -infinity) (coerce -0.0 (type-contagion (abs x) -infinity))))
+          (is (eql (/ (- (abs x)) -infinity) (coerce +0.0 (type-contagion (- (abs x)) -infinity))))
+          (is (eql (/ (abs x) +infinity) (coerce +0.0 (type-contagion (abs x) +infinity))))
+          (is (eql (/ (- (abs x)) +infinity) (coerce -0.0 (type-contagion (- (abs x)) +infinity))))
+
+          (type-and-value-check (/ -infinity (abs x)) -infinity)
+          (type-and-value-check (/ -infinity (- (abs x))) +infinity)
+          (type-and-value-check (/ +infinity (abs x)) +infinity)
+          (type-and-value-check (/ +infinity (- (abs x))) -infinity))))
+
+(test ieee-fp.0021.infinity-1+/1-
+      (for-all-infinities -infinity +infinity
+        (type-and-value-check (1- -infinity) -infinity)
+        (type-and-value-check (1- +infinity) +infinity)
+        (type-and-value-check (1+ -infinity) -infinity)
+        (type-and-value-check (1+ +infinity) +infinity)))
+
+(test ieee-fp.0022.infinity-abs
+      (for-all-infinities -infinity +infinity
+        (type-and-value-check (abs -infinity) +infinity)
+        (type-and-value-check (abs +infinity) +infinity)))
+
+(test ieee-fp.0023.infinity-exp
+      (for-all-infinities -infinity +infinity
+        (type-and-value-check (exp -infinity) 0)
+        (type-and-value-check (exp +infinity) +infinity)))
+
+(test ieee-fp.0024.infinity-expt
+      (without-fpe-traps
+        (flet ((relaxed-eql (x y)
+                 ;; (relaxed-eql +0.0 -0.0) => true
+                 (and (eql (type-of x) (type-of y))
+                      (= x y))))
+          (for-all-infinities -infinity +infinity
+            ;; x^∞: |x| < 1
+            (for-all-number-subtypes (x (and real (not integer)) 1/2)
+              (type-and-value-check (expt x -infinity) +infinity)
+              #| (complex-equality #'eql (expt (- x) -infinity)
+              (coerce +infinity (type-contagion x -infinity)) ;
+              (coerce (ext:nan) (type-contagion x -infinity))) |#
+              (type-and-value-check (expt x +infinity) 0)
+              #| (complex-equality #'relaxed-eql (expt (- x) +infinity)
+              (coerce 0.0 (type-contagion x +infinity)) ;
+              (coerce 0.0 (type-contagion x +infinity))) |#)
+            ;; x^∞: |x| = 1
+            (for-all-number-subtypes (x (or fixnum float) 1)
+              (type-and-value-check (expt x -infinity) 1)
+              (without-fpe-traps
+                (complex-equality #'eql (expt (- x) -infinity)
+                                  (coerce (ext:nan) (type-contagion x -infinity))
+                                  (coerce (ext:nan) (type-contagion x -infinity))))
+              (type-and-value-check (expt x +infinity) 1)
+              (without-fpe-traps
+                (complex-equality #'eql (expt (- x) +infinity)
+                                  (coerce (ext:nan) (type-contagion x +infinity))
+                                  (coerce (ext:nan) (type-contagion x +infinity)))))
+            ;; x^∞: |x| > 1
+            (for-all-number-subtypes (x (or fixnum float) 2)
+              (type-and-value-check (expt x -infinity) 0)
+              #| (complex-equality #'relaxed-eql (expt (- x) -infinity)
+              (coerce 0.0 (type-contagion x -infinity)) ;
+              (coerce 0.0 (type-contagion x -infinity))) |#
+              (type-and-value-check (expt x +infinity) +infinity)
+              #| (complex-equality #'eql (expt (- x) +infinity)
+              (coerce +infinity (type-contagion x +infinity)) ;
+              (coerce (ext:nan) (type-contagion x +infinity))) |#)
+            ;; ∞^y: even y
+            (for-all-number-subtypes (y float 2)
+              (is (eql (realpart (expt -infinity y))
+                       (coerce +infinity (type-contagion -infinity y))))
+              (complex-equality #'relaxed-eql (expt -infinity (- y))
+                                (coerce 0.0 (type-contagion -infinity y))
+                                (coerce 0.0 (type-contagion -infinity y)))
+              (type-and-value-check (expt +infinity y) +infinity)
+              (type-and-value-check (expt +infinity (- y)) 0.0))
+            (for-all-number-subtypes (y fixnum 2)
+              (type-and-value-check (expt -infinity y) +infinity)
+              (type-and-value-check (expt -infinity (- y)) 0.0)
+              (type-and-value-check (expt +infinity y) +infinity)
+              (type-and-value-check (expt +infinity (- y)) 0.0))
+            ;; ∞^y: odd y
+            (for-all-number-subtypes (y float 3)
+              ;; Should be -infinity, but C99 complex functions return +infinity 
+              #| (is (eql (realpart (expt -infinity y))
+              (coerce -infinity (type-contagion -infinity y)))) |#
+              (complex-equality #'relaxed-eql (expt -infinity (- y))
+                                (coerce 0.0 (type-contagion -infinity y))
+                                (coerce 0.0 (type-contagion -infinity y)))
+              (type-and-value-check (expt +infinity y) +infinity)
+              (type-and-value-check (expt +infinity (- y)) 0.0))
+            (for-all-number-subtypes (y fixnum 3)
+              (type-and-value-check (expt -infinity y) -infinity)
+              (complex-equality #'relaxed-eql (expt -infinity (- y))
+                                (coerce 0.0 (type-contagion -infinity y))
+                                (coerce 0.0 (type-contagion -infinity y)))
+              (type-and-value-check (expt +infinity y) +infinity)
+              (type-and-value-check (expt +infinity (- y)) 0.0)) 
+            ;; ∞^y: non-integer y
+            (for-all-number-subtypes (y (and real (not integer)) 5/2)
+              (complex-equality #'relaxed-eql (expt -infinity (- y))
+                                (coerce 0.0 (type-contagion -infinity y))
+                                (coerce 0.0 (type-contagion -infinity y)))
+              (type-and-value-check (expt +infinity y) +infinity)
+              (type-and-value-check (expt +infinity (- y)) 0.0))
+            ;; ∞^∞
+            #| (complex-equality #'eql (expt -infinity -infinity)
+            (coerce 0.0 (type-contagion -infinity -infinity)) ;
+            (coerce -0.0 (type-contagion -infinity -infinity))) |#
+            #| (complex-equality #'eql (expt -infinity +infinity)
+            (coerce +infinity (type-contagion -infinity +infinity)) ;
+            (coerce (ext:nan) (type-contagion -infinity +infinity))) |#
+            (type-and-value-check (expt +infinity -infinity) 0.0)
+            (type-and-value-check (expt +infinity +infinity) +infinity)))))
+
+(test ieee-fp.0025.infinity-log
+      (for-all-infinities -infinity +infinity
+        (complex-equality #'approx= (log -infinity) +infinity pi)
+        (type-and-value-check (log +infinity) +infinity)))
+
+(test ieee-fp.0026.infinity-signum
+      (for-all-infinities -infinity +infinity
+        (type-and-value-check (signum -infinity) -1)
+        (type-and-value-check (signum +infinity) +1)))
+
+(test ieee-fp.0027.infinity-sqrt
+      (for-all-infinities -infinity +infinity
+        (without-fpe-traps
+          (is (= (imagpart (sqrt -infinity)) +infinity)))
+        (type-and-value-check (sqrt +infinity) +infinity)))
+
+(test ieee-fp.0028.infinity-cis
+      (without-fpe-traps
+        (for-all-infinities -infinity +infinity
+          (is (eql (cis -infinity)
+                   (complex (coerce (ext:nan) (type-of -infinity))
+                            (coerce (ext:nan) (type-of -infinity)))))
+          (is (eql (cis +infinity)
+                   (complex (coerce (ext:nan) (type-of +infinity))
+                            (coerce (ext:nan) (type-of +infinity))))))))
+
+(test ieee-fp.0029.infinity-conjugate
+      (for-all-infinities -infinity +infinity
+        (type-and-value-check (conjugate -infinity) -infinity)
+        (type-and-value-check (conjugate +infinity) +infinity)
+        #| (for-all-number-subtypes (x real)
+          (type-and-value-check (conjugate (complex x -infinity)) (complex x +infinity))
+          (type-and-value-check (conjugate (complex x +infinity)) (complex x -infinity))) |#))
+
+(test ieee-fp.0030.infinity-phase
+      (for-all-infinities -infinity +infinity
+        (is (approx= (phase -infinity) (- pi)))
+        (is (approx= (phase +infinity) 0))
+        #| (is (approx= (phase (complex 0 -infinity)) (- (/ pi 2))))
+        (is (approx= (phase (complex 0 +infinity)) (+ (/ pi 2))))
+        (is (approx= (phase (complex -infinity -infinity)) (- (* pi 3/4))))
+        (is (approx= (phase (complex -infinity +infinity)) (+ (* pi 3/4))))
+        (is (approx= (phase (complex +infinity -infinity)) (- (* pi 1/4))))
+        (is (approx= (phase (complex +infinity +infinity)) (+ (* pi 1/4)))) |#))
diff --git a/src/tests/universe.lisp b/src/tests/universe.lisp
index 112d81e5e..ff8184684 100644
--- a/src/tests/universe.lisp
+++ b/src/tests/universe.lisp
@@ -129,6 +129,12 @@
         ext:double-float-positive-infinity
         ext:long-float-positive-infinity))
 
+#+ieee-floating-point
+(defparameter *ieee-fp.nan*
+  (list (coerce (ext:nan) 'single-float)
+        (coerce (ext:nan) 'double-float)
+        (coerce (ext:nan) 'long-float)))
+
 (defparameter *ratios*
     '(1/3 1/1000 1/1000000000000000 -10/3 -1000/7 -987129387912381/13612986912361
       189729874978126783786123/1234678123487612347896123467851234671234))

From e385abb213a4b48ce213df384df8e5c51d914ef4 Mon Sep 17 00:00:00 2001
From: Marius Gerbershagen <marius.gerbershagen@gmail.com>
Date: Wed, 14 Aug 2019 21:46:39 +0200
Subject: [PATCH 10/10] ieee-fp: document the behaviour of numeric functions
 for infinity/NaN

---
 src/doc/manual/standards/numbers.txi | 43 ++++++++++++++++++++++++++--
 1 file changed, 41 insertions(+), 2 deletions(-)

diff --git a/src/doc/manual/standards/numbers.txi b/src/doc/manual/standards/numbers.txi
index d6f6eba93..049f5582e 100644
--- a/src/doc/manual/standards/numbers.txi
+++ b/src/doc/manual/standards/numbers.txi
@@ -5,6 +5,7 @@
 * Numbers - Numeric types::
 * Numbers - Floating point exceptions::
 * Numbers - Random-States::
+* Numbers - Infinity and Not a Number::
 * Numbers - Dictionary::
 * Numbers - C Reference::
 @end menu
@@ -85,8 +86,8 @@ a generalized boolean
 If @var{condition} is @code{last}, @var{flag} is ignored and the
 currently enabled floating point exceptions are returned in an
 implementation depended format (currently an integer). Otherwise,
-@var{flag} determines whether current thread will signal a floating
-point exception for the conditions passed in @var{condition}.
+@var{flag} determines whether the current thread will signal a
+floating point exception for the conditions passed in @var{condition}.
 @var{condition} can be either a symbol denoting a single condition,
 @code{t} for all conditions that are enabled by default or a value
 obtained from an earlier call to @code{ext:trap-fpe} with @code{last}.
@@ -106,6 +107,44 @@ random seed (an integer), an array of random bytes (mainly used for
 reading back printed random-state) and another random-state (syntactic
 sugar for copying the random-state).
 
+@node Numbers - Infinity and Not a Number
+@subsection Infinity and Not a Number
+The @ansi{} standard does not specify the behaviour of numeric
+functions for infinite or not number valued floating point numbers. If
+ECL is configured to support these special values (see the
+@code{--with-ieee-fp} configure option) and floating point exceptions
+are disabled, numeric functions generally return the same value as the
+corresponding C function. This means, that the output will be a NaN
+for a NaN input, and the ``mathematically correct'' value (which may
+be NaN, e.g. for ∞-∞) for an infinite real input. For complex floats,
+however, the return value of a numeric function called with a complex
+number for which the real or imaginary part is infinite, is
+undefined@footnote{The main reason for this is that some numeric
+functions for C complex numbers return mathematically incorrect
+values, for example sinh(i*∞) returns i*NaN instead of the
+mathematically correct i*∞. Keeping this consistent with our own
+implementation of complex arithmetic that is used when C complex
+numbers are not available would require to much work. Furthermore,
+complex arithmetic with infinities is unreliable anyway, since it
+quickly leads to NaN values (consider i*∞ = (0+i*1)*(∞+i*0) = NaN+i*∞;
+even this simple example is already mathematically incorrect).}.
+
+For other functions dealing with numbers, we adopt the following
+behaviour:
+
+@table @asis
+@item Comparison functions
+All numeric comparisons with @code{=,<,<=,>,>=} involving NaN return
+false. Comparing two NaNs of the same type with @code{eql} returns
+true.
+@item @code{min/max}
+NaN values are ignored, i.e. the maximum/minimum is taken only over
+the number valued parameters.
+@item Rounding functions
+All rounding functions signal an arithmetic-error if any of the given
+parameters are not number valued or infinite.
+@end table
+
 @node Numbers - Dictionary
 @subsection Dictionary