Merge pull request #531 from pydata/add_floordiv

Added floor division and log2

Author: FrancescAlted

ADDFUNCS.rst

@@ -54,7 +54,7 @@ Example:
     #define FUNC_FF(...)
     #endif
     ...
-    FUNC_FF(FUNC_MYFUNC_FF, "myfunc_ff", myfuncf, myfuncf2, vfMyfunc)
+    FUNC_FF(FUNC_MYFUNC_FF, "myfunc_ff", myfuncf, myfuncf2, vsMyfunc)
     FUNC_FF(FUNC_FF_LAST,    NULL,       NULL,    NULL,     NULL)
     #ifdef ELIDE_FUNC_FF
     #undef ELIDE_FUNC_FF
@@ -171,7 +171,7 @@ Add clauses to generate the FUNC_CODES from the ``functions.hpp`` header, making
     };
     #endif
 
-Some functions (e.g. ``fmod``, ``isnan``) are not available in MKL, and so must be hard-coded here as well:
+Some functions (e.g. ``fmod``, ``isnan``) are not available in MKL, and so must be hard-coded in ``bespoke_functions.hpp`` as well:
 
 .. code-block:: cpp
 
@@ -186,7 +186,7 @@ Some functions (e.g. ``fmod``, ``isnan``) are not available in MKL, and so must
     };
     #endif
 
-The complex case is slightlñy different (see other examples in the same file).
+The complex case is slightly different (see other examples in the same file).
 
 Add case handling to the ``check_program`` function
 

doc/user_guide.rst

@@ -188,10 +188,10 @@ Supported operators
 
 *NumExpr* supports the set of operators listed below:
 
-    * Bitwise operators (and, or, not, xor): :code:`&, |, ~, ^`
+    * Bitwise and logical operators (and, or, not, xor): :code:`&, |, ~, ^`
     * Comparison operators: :code:`<, <=, ==, !=, >=, >`
     * Unary arithmetic operators: :code:`-`
-    * Binary arithmetic operators: :code:`+, -, *, /, **, %, <<, >>`
+    * Binary arithmetic operators: :code:`+, -, *, /, //, **, %, <<, >>`
 
 
 Supported functions
@@ -203,22 +203,33 @@ The next are the current supported set:
       is true, number2 otherwise.
     * :code:`{isinf, isnan, isfinite}(float|complex): bool` -- returns element-wise True
       for ``inf`` or ``NaN``, ``NaN``, not ``inf`` respectively.
+    * :code:`signbit(float|complex): bool` -- returns element-wise True if signbit is set
+      False otherwise.
     * :code:`{sin,cos,tan}(float|complex): float|complex` -- trigonometric sine,
       cosine or tangent.
     * :code:`{arcsin,arccos,arctan}(float|complex): float|complex` -- trigonometric
       inverse sine, cosine or tangent.
     * :code:`arctan2(float1, float2): float` -- trigonometric inverse tangent of
       float1/float2.
+    * :code:`hypot(float1, float2): float` -- Euclidean distance between float1, float2
+    * :code:`nextafter(float1, float2): float` -- next representable floating-point value after
+      float1 in direction of float2
+    * :code:`copysign(float1, float2): float` -- return number with magnitude of float1 and
+      sign of float2
+    * :code:`{maximum,minimum}(float1, float2): float` -- return max/min of float1, float2
     * :code:`{sinh,cosh,tanh}(float|complex): float|complex` -- hyperbolic sine,
       cosine or tangent.
     * :code:`{arcsinh,arccosh,arctanh}(float|complex): float|complex` -- hyperbolic
       inverse sine, cosine or tangent.
-    * :code:`{log,log10,log1p}(float|complex): float|complex` -- natural, base-10 and
+    * :code:`{log,log10,log1p,log2}(float|complex): float|complex` -- natural, base-10 and
       log(1+x) logarithms.
     * :code:`{exp,expm1}(float|complex): float|complex` -- exponential and exponential
       minus one.
     * :code:`sqrt(float|complex): float|complex` -- square root.
-    * :code:`abs(float|complex): float|complex`  -- absolute value.
+    * :code:`trunc(float): float` -- round towards zero
+    * :code:`round(float|complex|int): float|complex|int` -- round to nearest integer (`rint`)
+    * :code:`sign(float|complex|int): float|complex|int` -- return -1, 0, +1 depending on sign
+    * :code:`abs(float|complex|int): float|complex|int`  -- absolute value.
     * :code:`conj(complex): complex` -- conjugate value.
     * :code:`{real,imag}(complex): float` -- real or imaginary part of complex.
     * :code:`complex(float, float): complex` -- complex from real and imaginary

numexpr/bespoke_functions.hpp

@@ -0,0 +1,292 @@
+#include <numpy/npy_cpu.h>
+#include <math.h>
+#include <string.h>
+#include <assert.h>
+#include <vector>
+#include "numexpr_config.hpp" // isnan definitions
+
+// Generic sign function
+inline int signi(int x) {return (0 < x) - (x < 0);}
+inline long signl(long x) {return (0 < x) - (x < 0);}
+inline double sign(double x){
+        // Floats: -1.0, 0.0, +1.0, NaN stays NaN
+        if (isnand(x)) {return NAN;}
+        if (x > 0) {return 1;}
+        if (x < 0) {return -1;}
+        return 0; // handles +0.0 and -0.0
+    }
+inline float signf(float x){
+        // Floats: -1.0, 0.0, +1.0, NaN stays NaN
+        if (isnanf_(x)) {return NAN;}
+        if (x > 0) {return 1;}
+        if (x < 0) {return -1;}
+        return 0; // handles +0.0 and -0.0
+    }
+
+// round function for ints
+inline int rinti(int x) {return x;}
+inline long rintl(long x) {return x;}
+// abs function for ints
+inline int fabsi(int x) {return x<0 ? -x: x;}
+inline long fabsl(long x) {return x<0 ? -x: x;}
+// fmod function for ints
+inline int fmodi(int x, int y) {return (int)fmodf((float)x, (float)y);}
+inline long fmodl(long x, long y)  {return (long)fmodf((long)x, (long)y);}
+
+#ifdef USE_VML
+static void viRint(MKL_INT n, const int* x, int* dest)
+{
+    memcpy(dest, x, n * sizeof(int)); // just copy x1 which is already int
+};
+
+static void vlRint(MKL_INT n, const long* x, long* dest)
+{
+    memcpy(dest, x, n * sizeof(long)); // just copy x1 which is already int
+};
+
+static void viFabs(MKL_INT n, const int* x, int* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j] = x[j] < 0 ? -x[j]: x[j];
+    };
+};
+
+static void vlFabs(MKL_INT n, const long* x, long* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j] = x[j] < 0 ? -x[j]: x[j];
+    };
+};
+
+/* Fake vsConj function just for casting purposes inside numexpr */
+static void vsConj(MKL_INT n, const float* x1, float* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j] = x1[j];
+    };
+};
+
+/* fmod not available in VML */
+static void vsfmod(MKL_INT n, const float* x1, const float* x2, float* dest)
+{
+    MKL_INT j;
+    for(j=0; j < n; j++) {
+    dest[j] = fmodf(x1[j], x2[j]);
+    };
+}
+static void vdfmod(MKL_INT n, const double* x1, const double* x2, double* dest)
+{
+    MKL_INT j;
+    for(j=0; j < n; j++) {
+    dest[j] = fmod(x1[j], x2[j]);
+    };
+};
+static void vifmod(MKL_INT n, const int* x1, const int* x2, int* dest)
+{
+    MKL_INT j;
+    for(j=0; j < n; j++) {
+    dest[j] = fmodi(x1[j], x2[j]);
+    };
+};
+static void vlfmod(MKL_INT n, const long* x1, const long* x2, long* dest)
+{
+    MKL_INT j;
+    for(j=0; j < n; j++) {
+    dest[j] = fmodl(x1[j], x2[j]);
+    };
+};
+
+/* no isnan, isfinite, isinf or signbit in VML */
+static void vsIsfinite(MKL_INT n, const float* x1, bool* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j] = isfinitef_(x1[j]);
+    };
+};
+static void vsIsinf(MKL_INT n, const float* x1, bool* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j] = isinff_(x1[j]);
+    };
+};
+static void vsIsnan(MKL_INT n, const float* x1, bool* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j] = isnanf_(x1[j]);
+    };
+};
+static void vsSignBit(MKL_INT n, const float* x1, bool* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j] = signbitf(x1[j]);
+    };
+};
+
+/* no isnan, isfinite, isinf, signbit in VML */
+static void vdIsfinite(MKL_INT n, const double* x1, bool* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j] = isfinited(x1[j]);
+    };
+};
+static void vdIsinf(MKL_INT n, const double* x1, bool* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j] = isinfd(x1[j]);
+    };
+};
+static void vdIsnan(MKL_INT n, const double* x1, bool* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j] = isnand(x1[j]);
+    };
+};
+static void vdSignBit(MKL_INT n, const double* x1, bool* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j] = signbit(x1[j]);
+    };
+};
+
+/* no isnan, isfinite or isinf in VML */
+static void vzIsfinite(MKL_INT n, const MKL_Complex16* x1, bool* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j] = isfinited(x1[j].real) && isfinited(x1[j].imag);
+    };
+};
+static void vzIsinf(MKL_INT n, const MKL_Complex16* x1, bool* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j] = isinfd(x1[j].real) || isinfd(x1[j].imag);
+    };
+};
+static void vzIsnan(MKL_INT n, const MKL_Complex16* x1, bool* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j] = isnand(x1[j].real) || isnand(x1[j].imag);
+    };
+};
+
+/* Fake vdConj function just for casting purposes inside numexpr */
+static void vdConj(MKL_INT n, const double* x1, double* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j] = x1[j];
+    };
+};
+
+/* various functions not available in VML */
+static void vzExpm1(MKL_INT n, const MKL_Complex16* x1, MKL_Complex16* dest)
+{
+    MKL_INT j;
+    vzExp(n, x1, dest);
+    for (j=0; j<n; j++) {
+        dest[j].real -= 1.0;
+    };
+};
+
+static void vzLog1p(MKL_INT n, const MKL_Complex16* x1, MKL_Complex16* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j].real = x1[j].real + 1;
+        dest[j].imag = x1[j].imag;
+    };
+    vzLn(n, dest, dest);
+};
+
+static void vzLog2(MKL_INT n, const MKL_Complex16* x1, MKL_Complex16* dest)
+{
+    MKL_INT j;
+    vzLn(n, x1, dest);
+    for (j=0; j<n; j++) {
+        dest[j].real = dest[j].real * M_LOG2_E;
+        dest[j].imag = dest[j].imag * M_LOG2_E;
+    };
+};
+
+static void vzRint(MKL_INT n, const MKL_Complex16* x1, MKL_Complex16* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j].real = rint(x1[j].real);
+        dest[j].imag = rint(x1[j].imag);
+    };
+};
+
+/* Use this instead of native vzAbs in VML as it seems to work badly */
+static void vzAbs_(MKL_INT n, const MKL_Complex16* x1, MKL_Complex16* dest)
+{
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        dest[j].real = sqrt(x1[j].real*x1[j].real + x1[j].imag*x1[j].imag);
+        dest[j].imag = 0;
+    };
+};
+
+/*sign functions*/
+static void vsSign(MKL_INT n, const float* x1, float* dest)
+{
+    MKL_INT j;
+    for(j=0; j < n; j++) {
+        dest[j] = signf(x1[j]);
+    };
+};
+static void vdSign(MKL_INT n, const double* x1, double* dest)
+{
+    MKL_INT j;
+    for(j=0; j < n; j++) {
+        dest[j] = sign(x1[j]);
+    };
+};
+static void viSign(MKL_INT n, const int* x1, int* dest)
+{
+    MKL_INT j;
+    for(j=0; j < n; j++) {
+        dest[j] = signi(x1[j]);
+    };
+};
+static void vlSign(MKL_INT n, const long* x1, long* dest)
+{
+    MKL_INT j;
+    for(j=0; j < n; j++) {
+        dest[j] = signl(x1[j]);
+    };
+};
+static void vzSign(MKL_INT n, const MKL_Complex16* x1, MKL_Complex16* dest)
+{
+    MKL_INT j;
+    double mag;
+    for(j=0; j < n; j++) {
+        mag = sqrt(x1[j].real*x1[j].real + x1[j].imag*x1[j].imag);
+        if (isnand(mag)) {
+            dest[j].real = NAN;
+            dest[j].imag = NAN;
+        }
+        else if (mag == 0) {
+            dest[j].real = 0;
+            dest[j].imag = 0;
+        }
+        else {
+            dest[j].real = x1[j].real / mag;
+            dest[j].imag = x1[j].imag / mag;
+        }
+    };
+};
+#endif