From 36a209cd5f648d4e7cdac062d8ecb321b561847d Mon Sep 17 00:00:00 2001
From: JeffBezanson <jeff.bezanson@gmail.com>
Date: Wed, 15 Apr 2009 23:54:43 +0000
Subject: [PATCH] making = a builtin fixing = and eqv? to work properly on NaNs
 fixing other comparison predicates to be consistent

---
 femtolisp/compiler.lsp |  5 +++--
 femtolisp/equal.c      | 20 +++++++++++------
 femtolisp/flisp.c      | 50 +++++++++++++++++++++++++++++++++++++++++-
 femtolisp/flisp.h      |  2 +-
 femtolisp/opcodes.h    |  2 +-
 femtolisp/system.lsp   |  7 +++---
 femtolisp/table.c      |  1 -
 femtolisp/unittest.lsp | 14 ++++++++++++
 llt/dtypes.h           |  1 +
 llt/hashing.c          |  4 ++--
 llt/operators.c        | 15 ++++++++-----
 llt/utils.h            |  3 ++-
 12 files changed, 99 insertions(+), 25 deletions(-)

diff --git a/femtolisp/compiler.lsp b/femtolisp/compiler.lsp
index 6f3dae9..7719f6f 100644
--- a/femtolisp/compiler.lsp
+++ b/femtolisp/compiler.lsp
@@ -17,7 +17,7 @@
     :cons :list :car :cdr :set-car! :set-cdr!
     :eval :apply
 
-    :+ :- :* :/ :< :compare
+    :+ :- :* :/ := :< :compare
 
     :vector :aref :aset! :for
 
@@ -40,7 +40,8 @@
 	 :set-cdr! 2      :eval     1
 	 :apply    2      :<        2
          :for      3      :compare  2
-         :aref     2      :aset!    3))
+         :aref     2      :aset!    3
+	 :=        2))
 
 (define 1/Instructions (table.invert Instructions))
 
diff --git a/femtolisp/equal.c b/femtolisp/equal.c
index 9f1f748..f7f40c8 100644
--- a/femtolisp/equal.c
+++ b/femtolisp/equal.c
@@ -34,17 +34,23 @@ static void eq_union(htable_t *table, value_t a, value_t b,
 }
 
 // a is a fixnum, b is a cprim
-static value_t compare_num_cprim(value_t a, value_t b, int eq)
+static value_t compare_num_cprim(value_t a, value_t b, int eq, int swap)
 {
     cprim_t *bcp = (cprim_t*)ptr(b);
     numerictype_t bt = cp_numtype(bcp);
     fixnum_t ia = numval(a);
     void *bptr = cp_data(bcp);
-    if (cmp_eq(&ia, T_FIXNUM, bptr, bt))
+    if (cmp_eq(&ia, T_FIXNUM, bptr, bt, 1))
         return fixnum(0);
     if (eq) return fixnum(1);
-    if (cmp_lt(&ia, T_FIXNUM, bptr, bt))
-        return fixnum(-1);
+    if (swap) {
+        if (cmp_lt(bptr, bt, &ia, T_FIXNUM))
+            return fixnum(-1);
+    }
+    else {
+        if (cmp_lt(&ia, T_FIXNUM, bptr, bt))
+            return fixnum(-1);
+    }
     return fixnum(1);
 }
 
@@ -87,7 +93,7 @@ static value_t bounded_compare(value_t a, value_t b, int bound, int eq)
             return (numval(a) < numval(b)) ? fixnum(-1) : fixnum(1);
         }
         if (iscprim(b)) {
-            return compare_num_cprim(a, b, eq);
+            return compare_num_cprim(a, b, eq, 0);
         }
         return fixnum(-1);
     case TAG_SYM:
@@ -104,7 +110,7 @@ static value_t bounded_compare(value_t a, value_t b, int bound, int eq)
             cprim_t *acp=(cprim_t*)ptr(a), *bcp=(cprim_t*)ptr(b);
             numerictype_t at=cp_numtype(acp), bt=cp_numtype(bcp);
             void *aptr=cp_data(acp), *bptr=cp_data(bcp);
-            if (cmp_eq(aptr, at, bptr, bt))
+            if (cmp_eq(aptr, at, bptr, bt, 1))
                 return fixnum(0);
             if (eq) return fixnum(1);
             if (cmp_lt(aptr, at, bptr, bt))
@@ -112,7 +118,7 @@ static value_t bounded_compare(value_t a, value_t b, int bound, int eq)
             return fixnum(1);
         }
         else if (isfixnum(b)) {
-            return fixnum(-numval(compare_num_cprim(b, a, eq)));
+            return compare_num_cprim(b, a, eq, 1);
         }
         break;
     case TAG_CVALUE:
diff --git a/femtolisp/flisp.c b/femtolisp/flisp.c
index 9b96b6b..473f497 100644
--- a/femtolisp/flisp.c
+++ b/femtolisp/flisp.c
@@ -68,7 +68,7 @@ static char *builtin_names[] =
       "eval", "apply",
 
       // arithmetic
-      "+", "-", "*", "/", "<", "compare",
+      "+", "-", "*", "/", "=", "<", "compare",
 
       // sequences
       "vector", "aref", "aset!", "for",
@@ -649,6 +649,33 @@ int isnumber(value_t v)
     return (isfixnum(v) || iscprim(v));
 }
 
+static int numeric_equals(value_t a, value_t b)
+{
+    value_t tmp;
+    if (isfixnum(b)) {
+        tmp=a; a=b; b=tmp;
+    }
+    void *aptr, *bptr;
+    numerictype_t at, bt;
+    if (!iscprim(b)) type_error("=", "number", b);
+    cprim_t *cp = (cprim_t*)ptr(b);
+    fixnum_t fv;
+    bt = cp_numtype(cp);
+    bptr = cp_data(cp);
+    if (isfixnum(a)) {
+        fv = numval(a);
+        at = T_FIXNUM;
+        aptr = &fv;
+    }
+    else if (iscprim(a)) {
+        cp = (cprim_t*)ptr(a);
+        at = cp_numtype(cp);
+        aptr = cp_data(cp);
+    }
+    else type_error("=", "number", a);
+    return cmp_eq(aptr, at, bptr, bt, 0);
+}
+
 // read -----------------------------------------------------------------------
 
 #include "read.c"
@@ -1289,6 +1316,16 @@ static value_t eval_sexpr(value_t e, value_t *penv, int tail, uint32_t envsz)
             argcount("compare", nargs, 2);
             v = compare(Stack[SP-2], Stack[SP-1]);
             break;
+        case F_NUMEQ:
+            argcount("=", nargs, 2);
+            v = Stack[SP-2]; e = Stack[SP-1];
+            if (bothfixnums(v, e)) {
+                v = (v == e) ? FL_T : FL_F;
+            }
+            else {
+                v = numeric_equals(v, e) ? FL_T : FL_F;
+            }
+            break;
         case F_LT:
             argcount("<", nargs, 2);
             if (bothfixnums(Stack[SP-2], Stack[SP-1])) {
@@ -1857,6 +1894,17 @@ static value_t apply_cl(uint32_t nargs)
                 PUSH(v);
             }
             break;
+        case F_NUMEQ:
+            v = Stack[SP-2]; e = Stack[SP-1];
+            if (bothfixnums(v, e)) {
+                v = (v == e) ? FL_T : FL_F;
+            }
+            else {
+                v = numeric_equals(v, e) ? FL_T : FL_F;
+            }
+            POPN(1);
+            Stack[SP-1] = v;
+            break;
         case OP_LT:
             if (bothfixnums(Stack[SP-2], Stack[SP-1])) {
                 v = (numval(Stack[SP-2]) < numval(Stack[SP-1])) ? FL_T : FL_F;
diff --git a/femtolisp/flisp.h b/femtolisp/flisp.h
index 0688587..fb12d01 100644
--- a/femtolisp/flisp.h
+++ b/femtolisp/flisp.h
@@ -125,7 +125,7 @@ enum {
 
     F_CONS, F_LIST, F_CAR, F_CDR, F_SETCAR, F_SETCDR,
     F_EVAL, F_APPLY,
-    F_ADD, F_SUB, F_MUL, F_DIV, F_LT, F_COMPARE,
+    F_ADD, F_SUB, F_MUL, F_DIV, F_NUMEQ, F_LT, F_COMPARE,
 
     F_VECTOR, F_AREF, F_ASET, F_FOR,
     F_TRUE, F_FALSE, F_NIL,
diff --git a/femtolisp/opcodes.h b/femtolisp/opcodes.h
index 9c8a01f..bcd4015 100644
--- a/femtolisp/opcodes.h
+++ b/femtolisp/opcodes.h
@@ -12,7 +12,7 @@ enum {
     OP_CONS, OP_LIST, OP_CAR, OP_CDR, OP_SETCAR, OP_SETCDR,
     OP_EVAL, OP_APPLY,
 
-    OP_ADD, OP_SUB, OP_MUL, OP_DIV, OP_LT, OP_COMPARE,
+    OP_ADD, OP_SUB, OP_MUL, OP_DIV, OP_NUMEQ, OP_LT, OP_COMPARE,
 
     OP_VECTOR, OP_AREF, OP_ASET, OP_FOR,
 
diff --git a/femtolisp/system.lsp b/femtolisp/system.lsp
index 3cb9683..a61caef 100644
--- a/femtolisp/system.lsp
+++ b/femtolisp/system.lsp
@@ -86,11 +86,10 @@
 	((eqv?       (caar lst) item) (car lst))
 	(#t          (assv item (cdr lst)))))
 
-(define =   eqv?)
-(define (/= a b) (not (eqv? a b)))
+(define (/= a b) (not (= a b)))
 (define (>  a b) (< b a))
-(define (<= a b) (not (< b a)))
-(define (>= a b) (not (< a b)))
+(define (<= a b) (or (< a b) (= a b)))
+(define (>= a b) (or (< b a) (= a b)))
 (define (negative? x) (< x 0))
 (define (zero? x)     (= x 0))
 (define (positive? x) (> x 0))
diff --git a/femtolisp/table.c b/femtolisp/table.c
index d2b10a1..b366ad8 100644
--- a/femtolisp/table.c
+++ b/femtolisp/table.c
@@ -171,7 +171,6 @@ value_t fl_table_foldl(value_t *args, uint32_t nargs)
     htable_t *h = totable(args[2], "table.foldl");
     size_t i, n = h->size;
     void **table = h->table;
-    value_t c;
     for(i=0; i < n; i+=2) {
         if (table[i+1] != HT_NOTFOUND) {
             args[1] = applyn(3, args[0],
diff --git a/femtolisp/unittest.lsp b/femtolisp/unittest.lsp
index 25cd37d..adfd1b4 100644
--- a/femtolisp/unittest.lsp
+++ b/femtolisp/unittest.lsp
@@ -69,6 +69,20 @@
 
 (assert (equal (string 'sym #byte(65) #wchar(945) "blah") "symA\u03B1blah"))
 
+; NaNs
+(assert (equal? +nan.0 +nan.0))
+(assert (not (= +nan.0 +nan.0)))
+(assert (not (= +nan.0 -nan.0)))
+(assert (equal? (< +nan.0 3) (> 3 +nan.0)))
+(assert (equal? (< +nan.0 (double 3)) (> (double 3) +nan.0)))
+(assert (equal? (< +nan.0 3) (> (double 3) +nan.0)))
+(assert (equal? (< +nan.0 (double 3)) (> 3 +nan.0)))
+(assert (equal? (< +nan.0 3) (< +nan.0 (double 3))))
+(assert (equal? (> +nan.0 3) (> +nan.0 (double 3))))
+(assert (equal? (< 3 +nan.0) (> +nan.0 (double 3))))
+(assert (equal? (> 3 +nan.0) (> (double 3) +nan.0)))
+(assert (not (>= +nan.0 +nan.0)))
+
 ; this crashed once
 (for 1 10 (lambda (i) 0))
 
diff --git a/llt/dtypes.h b/llt/dtypes.h
index 9fc7686..6dd154d 100644
--- a/llt/dtypes.h
+++ b/llt/dtypes.h
@@ -116,6 +116,7 @@ typedef u_ptrint_t uptrint_t;
 #define LABS(n) (((n)^((n)>>(NBITS-1))) - ((n)>>(NBITS-1)))
 #define NBABS(n,nb) (((n)^((n)>>((nb)-1))) - ((n)>>((nb)-1)))
 #define DFINITE(d) (((*(int64_t*)&(d))&0x7ff0000000000000LL)!=0x7ff0000000000000LL)
+#define DNAN(d) (((*(int64_t*)&(d))&0x7ff8000000000000LL)==0x7ff8000000000000LL)
 
 extern double D_PNAN;
 extern double D_NNAN;
diff --git a/llt/hashing.c b/llt/hashing.c
index 5129664..21c69f3 100644
--- a/llt/hashing.c
+++ b/llt/hashing.c
@@ -93,11 +93,11 @@ void llt_init()
     ios_init_stdstreams();
 
     D_PNAN = strtod("+NaN",NULL);
-    D_NNAN = strtod("-NaN",NULL);
+    D_NNAN = -strtod("+NaN",NULL);
     D_PINF = strtod("+Inf",NULL);
     D_NINF = strtod("-Inf",NULL);
     F_PNAN = strtof("+NaN",NULL);
-    F_NNAN = strtof("-NaN",NULL);
+    F_NNAN = -strtof("+NaN",NULL);
     F_PINF = strtof("+Inf",NULL);
     F_NINF = strtof("-Inf",NULL);
 }
diff --git a/llt/operators.c b/llt/operators.c
index 0d7440b..8acf054 100644
--- a/llt/operators.c
+++ b/llt/operators.c
@@ -235,16 +235,21 @@ int cmp_lt(void *a, numerictype_t atag, void *b, numerictype_t btag)
     return 0;
 }
 
-int cmp_eq(void *a, numerictype_t atag, void *b, numerictype_t btag)
+int cmp_eq(void *a, numerictype_t atag, void *b, numerictype_t btag,
+           int equalnans)
 {
-    if (atag==btag)
+    if (atag==btag && !equalnans)
         return cmp_same_eq(a, b, atag);
 
     double da = conv_to_double(a, atag);
     double db = conv_to_double(b, btag);
 
-    if ((int)atag >= T_FLOAT && (int)btag >= T_FLOAT)
+    if ((int)atag >= T_FLOAT && (int)btag >= T_FLOAT) {
+        if (equalnans && DNAN(da)) {
+            return *(uint64_t*)&da == *(uint64_t*)&db;
+        }
         return (da == db);
+    }
 
     if (da != db)
         return 0;
@@ -339,8 +344,8 @@ void test_operators()
     assert(cmp_lt(&d, T_DOUBLE, &i64, T_INT64));
     assert(!cmp_lt(&i64, T_INT64, &d, T_DOUBLE));
 
-    assert(!cmp_eq(&d, T_DOUBLE, &i64, T_INT64));
+    assert(!cmp_eq(&d, T_DOUBLE, &i64, T_INT64, 0));
     i64 = DBL_MAXINT;
-    assert(cmp_eq(&d, T_DOUBLE, &i64, T_INT64));
+    assert(cmp_eq(&d, T_DOUBLE, &i64, T_INT64, 0));
 }
 #endif
diff --git a/llt/utils.h b/llt/utils.h
index 2d9d16f..bdece65 100644
--- a/llt/utils.h
+++ b/llt/utils.h
@@ -75,7 +75,8 @@ uint32_t conv_to_uint32(void *data, numerictype_t tag);
 int cmp_same_lt(void *a, void *b, numerictype_t tag);
 int cmp_same_eq(void *a, void *b, numerictype_t tag);
 int cmp_lt(void *a, numerictype_t atag, void *b, numerictype_t btag);
-int cmp_eq(void *a, numerictype_t atag, void *b, numerictype_t btag);
+int cmp_eq(void *a, numerictype_t atag, void *b, numerictype_t btag,
+           int equalnans);
 
 #ifdef ARCH_X86_64
 #  define LEGACY_REGS "=Q"