Clarify logic for nan comparisons.

Some coding style cleanups in arraytypes.c.src.
Add documentation to release notes for 1.4.0.
Add documentation of sort order to the sort docstring.

3 files changed, 116 insertions, 39 deletions

diff --git a/doc/release/1.4.0-notes.rst b/doc/release/1.4.0-notes.rst
index 6a15e318b..e2df0322d 100644
--- a/doc/release/1.4.0-notes.rst
+++ b/doc/release/1.4.0-notes.rst
@@ -59,6 +59,9 @@ Enhancements
 
        Complex numbers with the same nan placements are sorted according to
        the non-nan part if it exists.
+    #. The type comparison functions have been made consistent with the new
+       sort order of nans. Searchsorted now works with sorted arrays
+       containing nan values.
     #. Complex division has been made more resistent to overflow.
     #. Complex floor division has been made more resistent to overflow.
 
diff --git a/numpy/core/fromnumeric.py b/numpy/core/fromnumeric.py
index 176ef3e6f..f7f584d3d 100644
--- a/numpy/core/fromnumeric.py
+++ b/numpy/core/fromnumeric.py
@@ -449,6 +449,22 @@ def sort(a, axis=-1, kind='quicksort', order=None):
     the last axis is faster and uses less space than sorting along
     any other axis.
 
+    The sort order for complex numbers is lexicographic. If both the real
+    and imaginary parts are non-nan then the order is determined by the
+    real parts except when they are equal, in which case the order is
+    determined by the imaginary parts.
+
+    Previous to numpy 1.4.0 sorting real and complex arrays containing nan
+    values led to undefined behaviour. In numpy versions >= 1.4.0 nan
+    values are sorted to the end. The extended sort order is:
+   
+         Real: [R, nan]
+         Complex: [R + Rj, R + nanj, nan + Rj, nan + nanj]
+   
+    where R is a non-nan real value. Complex values with the same nan
+    placements are sorted according to the non-nan part if it exists.
+    Non-nan values are sorted as before.
+
     Examples
     --------
     >>> a = np.array([[1,4],[3,1]])
@@ -528,6 +544,9 @@ def argsort(a, axis=-1, kind='quicksort', order=None):
     -----
     See `sort` for notes on the different sorting algorithms.
 
+    As of Numpy 1.4.0 argsort works with real/complex arrays containing
+    nan values. The enhanced sort order is documented in the numpy.sort.
+
     Examples
     --------
     One dimensional array:
@@ -664,6 +683,9 @@ def searchsorted(a, v, side='left'):
     -----
     Binary search is used to find the required insertion points.
 
+    As of Numpy 1.4.0 searchsorted works with real/complex arrays containing
+    nan values. The enhanced sort order is documented in the numpy.sort.
+
     Examples
     --------
     >>> np.searchsorted([1,2,3,4,5], 3)
diff --git a/numpy/core/src/multiarray/arraytypes.c.src b/numpy/core/src/multiarray/arraytypes.c.src
index 88d04f099..58bad98e1 100644
--- a/numpy/core/src/multiarray/arraytypes.c.src
+++ b/numpy/core/src/multiarray/arraytypes.c.src
@@ -1747,7 +1747,22 @@ VOID_nonzero (char *ip, PyArrayObject *ap)
 #undef __ALIGNED
 
 
-/****************** compare **********************************/
+/*
+ *****************************************************************************
+ **                         COMPARISON FUNCTIONS                            **
+ *****************************************************************************
+ */
+
+/* boolean type */
+
+static int
+BOOL_compare(Bool *ip1, Bool *ip2, PyArrayObject *NPY_UNUSED(ap))
+{
+    return (*ip1 ? (*ip2 ? 0 : 1) : (*ip2 ? -1 : 0));
+}
+
+
+/* integer types */
 
 /**begin repeat
  * #TYPE = BYTE, UBYTE, SHORT, USHORT, INT, UINT, LONG, ULONG,
@@ -1768,11 +1783,31 @@ static int
 /**end repeat**/
 
 
+/* float types */
+
+/*
+ * The real/complex comparison functions are compatible with the new sort
+ * order for nans introduced in numpy 1.4.0. All nan values now compare
+ * larger than non-nan values and are sorted to the end. The comparison
+ * order is:
+ *
+ *      Real: [R, nan]
+ *      Complex: [R + Rj, R + nanj, nan + Rj, nan + nanj]
+ *
+ *  where complex values with the same nan placements are sorted according
+ *  to the non-nan part if it exists. If both the real and imaginary parts
+ *  of complex types are non-nan the order is the same as the real parts
+ *  unless they happen to be equal, in which case the order is that of the
+ *  imaginary parts.
+ */
+
 /**begin repeat
  * #TYPE = FLOAT, DOUBLE, LONGDOUBLE#
  * #type = float, double, longdouble#
  */
 
+#define LT(a,b) ((a) < (b) || ((b) != (b) && (a) ==(a)))
+
 static int
 @TYPE@_compare(@type@ *pa, @type@ *pb)
 {
@@ -1780,10 +1815,10 @@ static int
     const @type@ b = *pb;
     int ret;
 
-    if (a < b || (b != b && a == a)) {
+    if (LT(a,b)) {
         ret = -1;
     }
-    else if (a > b || (a != a && b == b)) {
+    else if (LT(b,a)) {
         ret = 1;
     }
     else {
@@ -1810,19 +1845,19 @@ C@TYPE@_compare(@type@ *pa, @type@ *pb)
             ret = 1;
         }
     }
-    else if (ar > br) {
-        if (bi != bi && ai == ai) {
-            ret = -1;
+    else if (br < ar) {
+        if (bi == bi || ai != ai) {
+            ret = 1;
         }
         else {
-            ret = 1;
+            ret = -1;
         }
     }
     else if (ar == br || (ar != ar && br != br)) {
-        if (ai < bi || (bi != bi && ai == ai)) {
+        if (LT(ai,bi)) {
             ret = -1;
         }
-        else if (ai > bi || (ai != ai && bi == bi)) {
+        else if (LT(bi,ai)) {
             ret = 1;
         }
         else {
@@ -1838,26 +1873,32 @@ C@TYPE@_compare(@type@ *pa, @type@ *pb)
 
     return ret;
 }
+
+#undef LT
+
 /**end repeat**/
 
 
-static int
-BOOL_compare(Bool *ip1, Bool *ip2, PyArrayObject *NPY_UNUSED(ap))
-{
-    return (*ip1 ? (*ip2 ? 0 : 1) : (*ip2 ? -1 : 0));
-}
+/* object type */
 
 static int
 OBJECT_compare(PyObject **ip1, PyObject **ip2, PyArrayObject *NPY_UNUSED(ap))
 {
     if ((*ip1 == NULL) || (*ip2 == NULL)) {
-        if (ip1 == ip2) return 1;
-        if (ip1 == NULL) return -1;
+        if (ip1 == ip2) {
+            return 1;
+        }
+        if (ip1 == NULL) {
+            return -1;
+        }
         return 1;
     }
     return PyObject_Compare(*ip1, *ip2);
 }
 
+
+/* string type */
+
 static int
 STRING_compare(char *ip1, char *ip2, PyArrayObject *ap)
 {
@@ -1874,33 +1915,38 @@ STRING_compare(char *ip1, char *ip2, PyArrayObject *ap)
     return 0;
 }
 
-/* taken from Python */
+
+/* unicode type */
+
 static int
 UNICODE_compare(PyArray_UCS4 *ip1, PyArray_UCS4 *ip2,
                 PyArrayObject *ap)
 {
-    int itemsize=ap->descr->elsize;
-    PyArray_UCS4 c1, c2;
-
-    if (itemsize < 0) return 0;
+    int itemsize = ap->descr->elsize;
 
+    if (itemsize < 0) {
+        return 0;
+    }
     while(itemsize-- > 0) {
-        c1 = *ip1++;
-        c2 = *ip2++;
-
-        if (c1 != c2)
+        PyArray_UCS4 c1 = *ip1++;
+        PyArray_UCS4 c2 = *ip2++;
+        if (c1 != c2) {
             return (c1 < c2) ? -1 : 1;
+        }
     }
     return 0;
 }
 
-/* If fields are defined, then compare on first field and if equal
-   compare on second field.  Continue until done or comparison results
-   in not_equal.
 
-   Must align data passed on to sub-comparisons.
-*/
+/* void type */
 
+/*
+ * If fields are defined, then compare on first field and if equal
+ * compare on second field.  Continue until done or comparison results
+ * in not_equal.
+ *
+ * Must align data passed on to sub-comparisons.
+ */
 static int
 VOID_compare(char *ip1, char *ip2, PyArrayObject *ap)
 {
@@ -1908,17 +1954,18 @@ VOID_compare(char *ip1, char *ip2, PyArrayObject *ap)
     PyObject *names, *key;
     PyObject *tup, *title;
     char *nip1, *nip2;
-    int i, offset, res=0;
+    int i, offset, res = 0;
 
-    if (!PyArray_HASFIELDS(ap))
+    if (!PyArray_HASFIELDS(ap)) {
         return STRING_compare(ip1, ip2, ap);
-
+    }
     descr = ap->descr;
-    /* Compare on the first-field.  If equal, then
-       compare on the second-field, etc.
+    /*
+     * Compare on the first-field.  If equal, then
+     * compare on the second-field, etc.
      */
     names = descr->names;
-    for (i=0; i<PyTuple_GET_SIZE(names); i++) {
+    for (i = 0; i < PyTuple_GET_SIZE(names); i++) {
         key = PyTuple_GET_ITEM(names, i);
         tup = PyDict_GetItem(descr->fields, key);
         if (!PyArg_ParseTuple(tup, "Oi|O", &new, &offset,
@@ -1932,15 +1979,18 @@ VOID_compare(char *ip1, char *ip2, PyArrayObject *ap)
             if (((intp)(nip1) % new->alignment) != 0) {
                 /* create buffer and copy */
                 nip1 = _pya_malloc(new->elsize);
-                if (nip1 == NULL) goto finish;
+                if (nip1 == NULL) {
+                    goto finish;
+                }
                 memcpy(nip1, ip1+offset, new->elsize);
             }
             if (((intp)(nip2) % new->alignment) != 0) {
                 /* copy data to a buffer */
                 nip2 = _pya_malloc(new->elsize);
                 if (nip2 == NULL) {
-                    if (nip1 != ip1+offset)
+                    if (nip1 != ip1+offset) {
                         _pya_free(nip1);
+                    }
                     goto finish;
                 }
                 memcpy(nip2, ip2+offset, new->elsize);
@@ -1955,7 +2005,9 @@ VOID_compare(char *ip1, char *ip2, PyArrayObject *ap)
                 _pya_free(nip2);
             }
         }
-        if (res != 0) break;
+        if (res != 0) {
+            break;
+        }
     }
 
 finish: