Add fmax, fmin functions for floats and complex. These conform to the ieee

standard for these functions. For numpy they are extended to complex numbers,
where a complex is considered a nan if either the real or complex part is. In
this case nan + nan*1j is returned if both arguments are nans.

Modify maximum and minimum so that if either input is nan, then nan is returned.
Likewise for the complex version.

2 files changed, 72 insertions, 14 deletions

diff --git a/numpy/core/code_generators/generate_umath.py b/numpy/core/code_generators/generate_umath.py
index d831ceda9..533c3a037 100644
--- a/numpy/core/code_generators/generate_umath.py
+++ b/numpy/core/code_generators/generate_umath.py
@@ -6,6 +6,7 @@ sys.path.pop(0)
 
 Zero = "PyUFunc_Zero"
 One = "PyUFunc_One"
+Nan = "PyUFunc_Nan"
 None_ = "PyUFunc_None"
 
 class TypeDescription(object):
@@ -316,6 +317,16 @@ defdict = {
           TD(noobj),
           TD(O, f='_npy_ObjectMin')
           ),
+'fmax' :
+    Ufunc(2, 1, Nan,
+          "",
+          TD(inexact),
+          ),
+'fmin' :
+    Ufunc(2, 1, Nan,
+          "",
+          TD(inexact),
+          ),
 'bitwise_and' :
     Ufunc(2, 1, One,
           docstrings.get('numpy.core.umath.bitwise_and'),
diff --git a/numpy/core/src/umathmodule.c.src b/numpy/core/src/umathmodule.c.src
index 6f41cd4fa..1af330a15 100644
--- a/numpy/core/src/umathmodule.c.src
+++ b/numpy/core/src/umathmodule.c.src
@@ -1579,7 +1579,23 @@ static void
 
 /**begin repeat1
  * #kind = maximum, minimum#
- * #OP =  >, <#
+ * #OP =  >=, <=#
+ **/
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+    /*  */
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        *((@type@ *)op) = (in1 @OP@ in2 || isnan(in1)) ? in1 : in2;
+    }
+}
+/**end repeat1**/
+
+/**begin repeat1
+ * #kind = fmax, fmin#
+ * #OP =  >=, <=#
  **/
 static void
 @TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
@@ -1588,7 +1604,7 @@ static void
     BINARY_LOOP {
         const @type@ in1 = *(@type@ *)ip1;
         const @type@ in2 = *(@type@ *)ip2;
-        *((@type@ *)op) = in1 @OP@ in2 ? in1 : in2;
+        *((@type@ *)op) = (in1 @OP@ in2 || isnan(in2)) ? in1 : in2;
     }
 }
 /**end repeat1**/
@@ -1680,15 +1696,7 @@ static void
     /*  */
     UNARY_LOOP {
         const @type@ in1 = *(@type@ *)ip1;
-        if (in1 > 0) {
-            *((@type@ *)op) = 1;
-        }
-        else if (in1 < 0) {
-            *((@type@ *)op) = -1;
-        }
-        else {
-            *((@type@ *)op) = 0;
-        }
+        *((@type@ *)op) = in1 > 0 ? 1 : (in1 < 0 ? -1 : 0);
     }
 }
 
@@ -1988,9 +1996,12 @@ static void
     }
 }
 
+#define GE(xr,xi,yr,yi) (xr > yr || (xr == yr && xi >= yi))
+#define LE(xr,xi,yr,yi) (xr < yr || (xr == yr && xi <= yi))
 /**begin repeat1
  * #kind = maximum, minimum#
- * #OP = >, <#
+ * #OP1 = GE, LE#
+ * #OP2 = LE, GE#
  */
 static void
 @CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
@@ -2000,10 +2011,44 @@ static void
         const @type@ in1i = ((@type@ *)ip1)[1];
         const @type@ in2r = ((@type@ *)ip2)[0];
         const @type@ in2i = ((@type@ *)ip2)[1];
-        if (in1r @OP@ in2r || ((in1r == in2r) && (in1i @OP@ in2i))) {
+        if (@OP1@(in1r, in1i, in2r, in2i)) {
             ((@type@ *)op)[0] = in1r;
             ((@type@ *)op)[1] = in1i;
         }
+        else if (@OP2@(in1r, in1i, in2r, in2i)) {
+            ((@type@ *)op)[0] = in2r;
+            ((@type@ *)op)[1] = in2i;
+        }
+        else {
+            ((@type@ *)op)[0] = NAN;
+            ((@type@ *)op)[1] = 0;
+        }
+    }
+}
+/**end repeat1**/
+
+/**begin repeat1
+ * #kind = fmax, fmin#
+ * #OP1 = GE, LE#
+ */
+static void
+@CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+    BINARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        const @type@ in2r = ((@type@ *)ip2)[0];
+        const @type@ in2i = ((@type@ *)ip2)[1];
+        if (@OP1@(in1r, in1i, in2r, in2i) || isnan(in2r) || isnan(in2i)) {
+            if (isnan(in1r) || isnan(in1i)) {
+                ((@type@ *)op)[0] = NAN;
+                ((@type@ *)op)[1] = 0;
+            }
+            else {
+                ((@type@ *)op)[0] = in1r;
+                ((@type@ *)op)[1] = in1i;
+            }
+        }
         else {
             ((@type@ *)op)[0] = in2r;
             ((@type@ *)op)[1] = in2i;
@@ -2011,6 +2056,8 @@ static void
     }
 }
 /**end repeat1**/
+#undef GE
+#undef LE
 
 #define @CTYPE@_true_divide @CTYPE@_divide
 /**end repeat**/
@@ -2251,7 +2298,7 @@ PyMODINIT_FUNC initumath(void) {
 
     pinf = pinf_init();
     pzero = pzero_init();
-    mynan = pinf / pinf;
+    mynan = copysign(pinf / pinf, 1);
 
     PyModule_AddObject(m, "PINF", PyFloat_FromDouble(pinf));
     PyModule_AddObject(m, "NINF", PyFloat_FromDouble(-pinf));