[Numpy-svn] r5894 - in branches/ufunc_cleanup/numpy/core: code_generators src

[numpy-svn at scipy.org]

Author: charris
Date: 2008-10-02 21:19:17 -0500 (Thu, 02 Oct 2008)
New Revision: 5894

Modified:
   branches/ufunc_cleanup/numpy/core/code_generators/generate_umath.py
   branches/ufunc_cleanup/numpy/core/src/umathmodule.c.src
Log:
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.


Modified: branches/ufunc_cleanup/numpy/core/code_generators/generate_umath.py
===================================================================
--- branches/ufunc_cleanup/numpy/core/code_generators/generate_umath.py	2008-10-03 02:11:50 UTC (rev 5893)
+++ branches/ufunc_cleanup/numpy/core/code_generators/generate_umath.py	2008-10-03 02:19:17 UTC (rev 5894)
@@ -6,6 +6,7 @@
 
 Zero = "PyUFunc_Zero"
 One = "PyUFunc_One"
+Nan = "PyUFunc_Nan"
 None_ = "PyUFunc_None"
 
 class TypeDescription(object):
@@ -316,6 +317,16 @@
           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'),

Modified: branches/ufunc_cleanup/numpy/core/src/umathmodule.c.src
===================================================================
--- branches/ufunc_cleanup/numpy/core/src/umathmodule.c.src	2008-10-03 02:11:50 UTC (rev 5893)
+++ branches/ufunc_cleanup/numpy/core/src/umathmodule.c.src	2008-10-03 02:19:17 UTC (rev 5894)
@@ -1579,7 +1579,7 @@
 
 /**begin repeat1
  * #kind = maximum, minimum#
- * #OP =  >, <#
+ * #OP =  >=, <=#
  **/
 static void
 @TYPE at _@kind@(char **args, intp *dimensions, intp *steps, void *func)
@@ -1588,12 +1588,28 @@
     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(in1)) ? in1 : in2;
     }
 }
 /**end repeat1**/
 
+/**begin repeat1
+ * #kind = fmax, fmin#
+ * #OP =  >=, <=#
+ **/
 static void
+ at TYPE@_ at 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(in2)) ? in1 : in2;
+    }
+}
+/**end repeat1**/
+
+static void
 @TYPE at _floor_divide(char **args, intp *dimensions, intp *steps, void *func)
 {
     BINARY_LOOP {
@@ -1680,15 +1696,7 @@
     /*  */
     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 @@
     }
 }
 
+#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 at _@kind@(char **args, intp *dimensions, intp *steps, void *func)
@@ -2000,17 +2011,53 @@
         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
+ at CTYPE@_ at 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;
         }
     }
 }
 /**end repeat1**/
+#undef GE
+#undef LE
 
 #define @CTYPE at _true_divide @CTYPE at _divide
 /**end repeat**/
@@ -2251,7 +2298,7 @@
 
     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));