chore: import upstream snapshot with attribution

numexpr/__init__.py

@@ -0,0 +1,66 @@
+###################################################################
+#  Numexpr - Fast numerical array expression evaluator for NumPy.
+#
+#      License: MIT
+#      Author:  See AUTHORS.txt
+#
+#  See LICENSE.txt and LICENSES/*.txt for details about copyright and
+#  rights to use.
+####################################################################
+
+"""
+Numexpr is a fast numerical expression evaluator for NumPy.  With it,
+expressions that operate on arrays (like "3*a+4*b") are accelerated
+and use less memory than doing the same calculation in Python.
+
+See:
+
+https://github.com/pydata/numexpr
+
+for more info about it.
+
+"""
+
+from numexpr.interpreter import __BLOCK_SIZE1__, MAX_THREADS, use_vml
+
+is_cpu_amd_intel = False # DEPRECATION WARNING: WILL BE REMOVED IN FUTURE RELEASE
+
+# cpuinfo imports were moved into the test submodule function that calls them
+# to improve import times.
+
+from numexpr.expressions import E
+from numexpr.necompiler import (NumExpr, disassemble, evaluate, re_evaluate,
+                                validate)
+from numexpr.utils import (_init_num_threads, detect_number_of_cores,
+                           detect_number_of_threads, get_num_threads,
+                           get_vml_version, set_num_threads,
+                           set_vml_accuracy_mode, set_vml_num_threads)
+
+# Detect the number of cores
+ncores = detect_number_of_cores()
+# Initialize the number of threads to be used
+nthreads = _init_num_threads()
+# The default for VML is 1 thread (see #39)
+# set_vml_num_threads(1)
+
+from . import version
+
+__version__ = version.version
+
+def print_versions():
+    """Print the versions of software that numexpr relies on."""
+    try:
+        import numexpr.tests
+        return numexpr.tests.print_versions()
+    except ImportError:
+        # To maintain Python 2.6 compatibility we have simple error handling
+        raise ImportError('`numexpr.tests` could not be imported, likely it was excluded from the distribution.')
+
+def test(verbosity=1):
+    """Run all the tests in the test suite."""
+    try:
+        import numexpr.tests
+        return numexpr.tests.test(verbosity=verbosity)
+    except ImportError:
+        # To maintain Python 2.6 compatibility we have simple error handling
+        raise ImportError('`numexpr.tests` could not be imported, likely it was excluded from the distribution.')

numexpr/bespoke_functions.hpp

@@ -0,0 +1,339 @@
+#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
+// TODO: Have to add FUNC_III, FUNC_LLL signatures to functions.hpp to enable these
+// 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
+// To match Numpy behaviour for NaNs
+static void vsFmax_(MKL_INT n, const float* x1, const float* x2, float* dest)
+{
+    vsFmax(n, x1, x2, dest);
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        if (isnanf_(x1[j]) | isnanf_(x2[j])){
+            dest[j] = NAN;
+        }
+    };
+};
+
+static void vsFmin_(MKL_INT n, const float* x1, const float* x2, float* dest)
+{
+    vsFmin(n, x1, x2, dest);
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        if (isnanf_(x1[j]) | isnanf_(x2[j])){
+            dest[j] = NAN;
+        }
+    };
+};
+// To match Numpy behaviour for NaNs
+static void vdFmax_(MKL_INT n, const double* x1, const double* x2, double* dest)
+{
+    vdFmax(n, x1, x2, dest);
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        if (isnand(x1[j]) | isnand(x2[j])){
+            dest[j] = NAN;
+        }
+    };
+};
+
+static void vdFmin_(MKL_INT n, const double* x1, const double* x2, double* dest)
+{
+    vdFmin(n, x1, x2, dest);
+    MKL_INT j;
+    for (j=0; j<n; j++) {
+        if (isnand(x1[j]) | isnand(x2[j])){
+            dest[j] = NAN;
+        }
+    };
+};
+
+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]);
+    };
+};
+// TODO: Have to add FUNC_III, FUNC_LLL signatures to functions.hpp
+// 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

numexpr/complex_functions.hpp

@@ -0,0 +1,498 @@
+#ifndef NUMEXPR_COMPLEX_FUNCTIONS_HPP
+#define NUMEXPR_COMPLEX_FUNCTIONS_HPP
+
+/*********************************************************************
+  Numexpr - Fast numerical array expression evaluator for NumPy.
+
+      License: MIT
+      Author:  See AUTHORS.txt
+
+  See LICENSE.txt for details about copyright and rights to use.
+**********************************************************************/
+
+// Replace npy_cdouble with std::complex<double>
+#include <math.h> // NAN
+#include <complex>
+
+/* constants */
+static std::complex<double> nc_1(1., 0.);
+static std::complex<double> nc_half(0.5, 0.);
+static std::complex<double> nc_i(0., 1.);
+static std::complex<double> nc_i2(0., 0.5);
+/*
+static std::complex<double> nc_mi = {0., -1.};
+static std::complex<double> nc_pi2 = {M_PI/2., 0.};
+*/
+
+/* *************************** WARNING *****************************
+Due to the way Numexpr places the results of operations, the *x and *r
+pointers do point to the same address (apparently this doesn't happen
+in NumPy).  So, measures should be taken so as to not to reuse *x
+after the first *r has been overwritten.
+*********************************************************************
+*/
+
+static void
+nc_assign(std::complex<double> *x, std::complex<double> *r)
+{
+  r->real(x->real());
+  r->imag(x->imag());
+  return;
+}
+
+static void
+nc_sum(std::complex<double> *a, std::complex<double> *b, std::complex<double> *r)
+{
+    r->real(a->real() + b->real());
+    r->imag(a->imag() + b->imag());
+    return;
+}
+
+static void
+nc_diff(std::complex<double> *a, std::complex<double> *b, std::complex<double> *r)
+{
+    r->real(a->real() - b->real());
+    r->imag(a->imag() - b->imag());
+    return;
+}
+
+static void
+nc_neg(std::complex<double> *a, std::complex<double> *r)
+{
+    r->real(-a->real());
+    r->imag(-a->imag());
+    return;
+}
+
+static void
+nc_conj(std::complex<double> *a, std::complex<double> *r)
+{
+    r->real(a->real());
+    r->imag(-a->imag());
+    return;
+}
+
+// Needed for allowing the internal casting in numexpr machinery for
+// conjugate operations
+inline float fconjf(float x)
+{
+    return x;
+}
+
+// Needed for allowing the internal casting in numexpr machinery for
+// conjugate operations
+inline double fconj(double x)
+{
+    return x;
+}
+
+static void
+nc_prod(std::complex<double> *a, std::complex<double> *b, std::complex<double> *r)
+{
+    double ar=a->real(), br=b->real(), ai=a->imag(), bi=b->imag();
+    r->real(ar*br - ai*bi);
+    r->imag(ar*bi + ai*br);
+    return;
+}
+
+static void
+nc_quot(std::complex<double> *a, std::complex<double> *b, std::complex<double> *r)
+{
+    double ar=a->real(), br=b->real(), ai=a->imag(), bi=b->imag();
+    double d = br*br + bi*bi;
+    r->real((ar*br + ai*bi)/d);
+    r->imag((ai*br - ar*bi)/d);
+    return;
+}
+
+static void
+nc_sqrt(std::complex<double> *x, std::complex<double> *r)
+{
+    double s,d;
+    if (x->real() == 0. && x->imag() == 0.)
+        *r = *x;
+    else {
+        s = sqrt((fabs(x->real()) + hypot(x->real(),x->imag()))/2);
+        d = x->imag()/(2*s);
+        if (x->real() > 0.) {
+            r->real(s);
+            r->imag(d);
+        }
+        else if (x->imag() >= 0.) {
+            r->real(d);
+            r->imag(s);
+        }
+        else {
+            r->real(-d);
+            r->imag(-s);
+        }
+    }
+    return;
+}
+
+static void
+nc_log(std::complex<double> *x, std::complex<double> *r)
+{
+    double l = hypot(x->real(),x->imag());
+    r->imag(atan2(x->imag(), x->real()));
+    r->real(log(l));
+    return;
+}
+
+static void
+nc_log1p(std::complex<double> *x, std::complex<double> *r)
+{
+    double l = hypot(x->real() + 1.0,x->imag());
+    r->imag(atan2(x->imag(), x->real() + 1.0));
+    r->real(log(l));
+    return;
+}
+
+static void
+nc_exp(std::complex<double> *x, std::complex<double> *r)
+{
+    double a = exp(x->real());
+    r->real(a*cos(x->imag()));
+    r->imag(a*sin(x->imag()));
+    return;
+}
+
+static void
+nc_expm1(std::complex<double> *x, std::complex<double> *r)
+{
+    double a = sin(x->imag() / 2);
+    double b = exp(x->real());
+    r->real(expm1(x->real()) * cos(x->imag()) - 2 * a * a);
+    r->imag(b * sin(x->imag()));
+    return;
+}
+
+static void
+nc_pow(std::complex<double> *a, std::complex<double> *b, std::complex<double> *r)
+{
+    npy_intp n;
+    double ar=a->real(), br=b->real(), ai=a->imag(), bi=b->imag();
+
+    if (br == 0. && bi == 0.) {
+        r->real(1.);
+        r->imag(0.);
+        return;
+    }
+    if (ar == 0. && ai == 0.) {
+        r->real(0.);
+        r->imag(0.);
+        return;
+    }
+    if (bi == 0 && (n=(npy_intp)br) == br) {
+        if (n > -100 && n < 100) {
+        std::complex<double> p, aa;
+        npy_intp mask = 1;
+        if (n < 0) n = -n;
+        aa = nc_1;
+        p.real(ar); p.imag(ai);
+        while (1) {
+            if (n & mask)
+                nc_prod(&aa,&p,&aa);
+            mask <<= 1;
+            if (n < mask || mask <= 0) break;
+            nc_prod(&p,&p,&p);
+        }
+        r->real(aa.real()); r->imag(aa.imag());
+        if (br < 0) nc_quot(&nc_1, r, r);
+        return;
+        }
+    }
+    /* complexobject.c uses an inline version of this formula
+       investigate whether this had better performance or accuracy */
+    nc_log(a, r);
+    nc_prod(r, b, r);
+    nc_exp(r, r);
+    return;
+}
+
+
+static void
+nc_prodi(std::complex<double> *x, std::complex<double> *r)
+{
+    double xr = x->real();
+    r->real(-x->imag());
+    r->imag(xr);
+    return;
+}
+
+
+static void
+nc_acos(std::complex<double> *x, std::complex<double> *r)
+{
+    std::complex<double> a, *pa=&a;
+
+    nc_assign(x, pa);
+    nc_prod(x,x,r);
+    nc_diff(&nc_1, r, r);
+    nc_sqrt(r, r);
+    nc_prodi(r, r);
+    nc_sum(pa, r, r);
+    nc_log(r, r);
+    nc_prodi(r, r);
+    nc_neg(r, r);
+    return;
+    /* return nc_neg(nc_prodi(nc_log(nc_sum(x,nc_prod(nc_i,
+       nc_sqrt(nc_diff(nc_1,nc_prod(x,x))))))));
+    */
+}
+
+static void
+nc_acosh(std::complex<double> *x, std::complex<double> *r)
+{
+    std::complex<double> t, a, *pa=&a;
+
+    nc_assign(x, pa);
+    nc_sum(x, &nc_1, &t);
+    nc_sqrt(&t, &t);
+    nc_diff(x, &nc_1, r);
+    nc_sqrt(r, r);
+    nc_prod(&t, r, r);
+    nc_sum(pa, r, r);
+    nc_log(r, r);
+    return;
+    /*
+      return nc_log(nc_sum(x,
+      nc_prod(nc_sqrt(nc_sum(x,nc_1)), nc_sqrt(nc_diff(x,nc_1)))));
+    */
+}
+
+static void
+nc_asin(std::complex<double> *x, std::complex<double> *r)
+{
+    std::complex<double> a, *pa=&a;
+    nc_prodi(x, pa);
+    nc_prod(x, x, r);
+    nc_diff(&nc_1, r, r);
+    nc_sqrt(r, r);
+    nc_sum(pa, r, r);
+    nc_log(r, r);
+    nc_prodi(r, r);
+    nc_neg(r, r);
+    return;
+    /*
+      return nc_neg(nc_prodi(nc_log(nc_sum(nc_prod(nc_i,x),
+      nc_sqrt(nc_diff(nc_1,nc_prod(x,x)))))));
+    */
+}
+
+
+static void
+nc_asinh(std::complex<double> *x, std::complex<double> *r)
+{
+    std::complex<double> a, *pa=&a;
+    nc_assign(x, pa);
+    nc_prod(x, x, r);
+    nc_sum(&nc_1, r, r);
+    nc_sqrt(r, r);
+    nc_sum(r, pa, r);
+    nc_log(r, r);
+    return;
+    /*
+      return nc_log(nc_sum(nc_sqrt(nc_sum(nc_1,nc_prod(x,x))),x));
+    */
+}
+
+static void
+nc_atan(std::complex<double> *x, std::complex<double> *r)
+{
+    std::complex<double> a, *pa=&a;
+    nc_diff(&nc_i, x, pa);
+    nc_sum(&nc_i, x, r);
+    nc_quot(r, pa, r);
+    nc_log(r,r);
+    nc_prod(&nc_i2, r, r);
+    return;
+    /*
+      return nc_prod(nc_i2,nc_log(nc_quot(nc_sum(nc_i,x),nc_diff(nc_i,x))));
+    */
+}
+
+static void
+nc_atanh(std::complex<double> *x, std::complex<double> *r)
+{
+    std::complex<double> a, b, *pa=&a, *pb=&b;
+    nc_assign(x, pa);
+    nc_diff(&nc_1, pa, r);
+    nc_sum(&nc_1, pa, pb);
+    nc_quot(pb, r, r);
+    nc_log(r, r);
+    nc_prod(&nc_half, r, r);
+    return;
+    /*
+      return nc_prod(nc_half,nc_log(nc_quot(nc_sum(nc_1,x),nc_diff(nc_1,x))));
+    */
+}
+
+static void
+nc_cos(std::complex<double> *x, std::complex<double> *r)
+{
+    double xr=x->real(), xi=x->imag();
+    r->real(cos(xr)*cosh(xi));
+    r->imag(-sin(xr)*sinh(xi));
+    return;
+}
+
+static void
+nc_cosh(std::complex<double> *x, std::complex<double> *r)
+{
+    double xr=x->real(), xi=x->imag();
+    r->real(cos(xi)*cosh(xr));
+    r->imag(sin(xi)*sinh(xr));
+    return;
+}
+
+
+#define M_LOG10_E 0.434294481903251827651128918916605082294397
+#define M_LOG2_E  1.44269504088896340735992468100189213742664
+
+
+static void
+nc_log10(std::complex<double> *x, std::complex<double> *r)
+{
+    nc_log(x, r);
+    r->real(r->real() * M_LOG10_E);
+    r->imag(r->imag() * M_LOG10_E);
+    return;
+}
+
+static void
+nc_log2(std::complex<double> *x, std::complex<double> *r)
+{
+    nc_log(x, r);
+    r->real(r->real() * M_LOG2_E);
+    r->imag(r->imag() * M_LOG2_E);
+    return;
+}
+
+static void
+nc_sin(std::complex<double> *x, std::complex<double> *r)
+{
+    double xr=x->real(), xi=x->imag();
+    r->real(sin(xr)*cosh(xi));
+    r->imag(cos(xr)*sinh(xi));
+    return;
+}
+
+static void
+nc_sinh(std::complex<double> *x, std::complex<double> *r)
+{
+    double xr=x->real(), xi=x->imag();
+    r->real(cos(xi)*sinh(xr));
+    r->imag(sin(xi)*cosh(xr));
+    return;
+}
+
+static void
+nc_tan(std::complex<double> *x, std::complex<double> *r)
+{
+    double xr = x->real();
+    double xi = x->imag();
+    double imag_part;
+
+    double denom = cos(2*xr) + cosh(2*xi);
+    // handle overflows
+    if (xi > 20) {
+        imag_part = 1.0 / (1.0 + exp(-4*xi));
+    } else if (xi < -20) {
+        imag_part = -1.0 / (1.0 + exp(4*xi));
+    } else {
+        imag_part = sinh(2*xi) / denom;
+    }
+    double real_part = sin(2*xr) / denom;
+
+    r->real(real_part);
+    r->imag(imag_part);
+    return;
+}
+
+static void
+nc_tanh(std::complex<double> *x, std::complex<double> *r)
+{
+    double xr = x->real();
+    double xi = x->imag();
+    double real_part;
+    double denom = cosh(2*xr) + cos(2*xi);
+    // handle overflows
+    if (xr > 20) {
+        real_part = 1.0 / (1.0 + exp(-4*xr));
+    } else if (xr < -20) {
+        real_part = -1.0 / (1.0 + exp(4*xr));
+    } else {
+        real_part = sinh(2*xr) / denom;
+    }
+    double imag_part = sin(2*xi) / denom;
+
+    r->real(real_part);
+    r->imag(imag_part);
+    return;
+}
+
+static void
+nc_abs(std::complex<double> *x, std::complex<double> *r)
+{
+    r->real(sqrt(x->real()*x->real() + x->imag()*x->imag()));
+    r->imag(0);
+}
+
+static void
+nc_rint(std::complex<double> *x, std::complex<double> *r)
+{
+    r->real(rint(x->real()));
+    r->imag(rint(x->imag()));
+}
+
+static bool
+nc_isinf(std::complex<double> *x)
+{
+    double xr=x->real(), xi=x->imag();
+    bool bi,br;
+    bi = isinfd(xi);
+    br = isinfd(xr);
+    return bi || br;
+}
+
+static bool
+nc_isnan(std::complex<double> *x)
+{
+    double xr=x->real(), xi=x->imag();
+    bool bi,br;
+    bi = isnand(xi);
+    br = isnand(xr);
+    return bi || br;
+}
+
+static bool
+nc_isfinite(std::complex<double> *x)
+{
+    double xr=x->real(), xi=x->imag();
+    bool bi,br;
+    bi = isfinited(xi);
+    br = isfinited(xr);
+    return bi && br;
+}
+
+static void
+nc_sign(std::complex<double> *x, std::complex<double> *r)
+{
+    if (nc_isnan(x)){
+        r->real(NAN);
+        r->imag(NAN);
+    }
+    std::complex<double> mag;
+    nc_abs(x, &mag);
+    if (mag.real() == 0){
+        r->real(0);
+        r->imag(0);
+    }
+    else{
+        r->real(x->real()/mag.real());
+        r->imag(x->imag()/mag.real());
+    }
+}
+
+#endif // NUMEXPR_COMPLEX_FUNCTIONS_HPP

numexpr/cpuinfo.py

@@ -0,0 +1,861 @@
+###################################################################
+#  cpuinfo - Get information about CPU
+#
+#      License: BSD
+#      Author:  Pearu Peterson <pearu@cens.ioc.ee>
+#
+#  See LICENSES/cpuinfo.txt for details about copyright and
+#  rights to use.
+####################################################################
+
+"""
+cpuinfo
+
+Copyright 2002 Pearu Peterson all rights reserved,
+Pearu Peterson <pearu@cens.ioc.ee>
+Permission to use, modify, and distribute this software is given under the
+terms of the NumPy (BSD style) license.  See LICENSE.txt that came with
+this distribution for specifics.
+
+NO WARRANTY IS EXPRESSED OR IMPLIED.  USE AT YOUR OWN RISK.
+Pearu Peterson
+"""
+
+__all__ = ['cpu']
+
+import inspect
+import os
+import platform
+import re
+import subprocess
+import sys
+import types
+import warnings
+
+is_cpu_amd_intel = False # DEPRECATION WARNING: WILL BE REMOVED IN FUTURE RELEASE
+
+def getoutput(cmd, successful_status=(0,), stacklevel=1):
+    try:
+        p = subprocess.Popen(cmd, stdout=subprocess.PIPE)
+        output, _ = p.communicate()
+        status = p.returncode
+    except EnvironmentError as e:
+        warnings.warn(str(e), UserWarning, stacklevel=stacklevel)
+        return False, ''
+    if os.WIFEXITED(status) and os.WEXITSTATUS(status) in successful_status:
+        return True, output
+    return False, output
+
+
+def command_info(successful_status=(0,), stacklevel=1, **kw):
+    info = {}
+    for key in kw:
+        ok, output = getoutput(kw[key], successful_status=successful_status,
+                               stacklevel=stacklevel + 1)
+        if ok:
+            info[key] = output.strip()
+    return info
+
+
+def command_by_line(cmd, successful_status=(0,), stacklevel=1):
+    ok, output = getoutput(cmd, successful_status=successful_status,
+                           stacklevel=stacklevel + 1)
+    if not ok:
+        return
+
+    # XXX: check
+    output = output.decode('ascii')
+
+    for line in output.splitlines():
+        yield line.strip()
+
+
+def key_value_from_command(cmd, sep, successful_status=(0,),
+                           stacklevel=1):
+    d = {}
+    for line in command_by_line(cmd, successful_status=successful_status,
+                                stacklevel=stacklevel + 1):
+        l = [s.strip() for s in line.split(sep, 1)]
+        if len(l) == 2:
+            d[l[0]] = l[1]
+    return d
+
+
+class CPUInfoBase(object):
+    """Holds CPU information and provides methods for requiring
+    the availability of various CPU features.
+    """
+
+    def _try_call(self, func):
+        try:
+            return func()
+        except:
+            pass
+
+    def __getattr__(self, name):
+        if not name.startswith('_'):
+            if hasattr(self, '_' + name):
+                attr = getattr(self, '_' + name)
+                if inspect.ismethod(attr):
+                    return lambda func=self._try_call, attr=attr: func(attr)
+            else:
+                return lambda: None
+        raise AttributeError(name)
+
+    def _getNCPUs(self):
+        return 1
+
+    def __get_nbits(self):
+        abits = platform.architecture()[0]
+        nbits = re.compile(r'(\d+)bit').search(abits).group(1)
+        return nbits
+
+    def _is_32bit(self):
+        return self.__get_nbits() == '32'
+
+    def _is_64bit(self):
+        return self.__get_nbits() == '64'
+
+
+class LinuxCPUInfo(CPUInfoBase):
+    info = None
+
+    def __init__(self):
+        if self.info is not None:
+            return
+        info = [{}]
+        ok, output = getoutput(['uname', '-m'])
+        if ok:
+            info[0]['uname_m'] = output.strip()
+        try:
+            fo = open('/proc/cpuinfo')
+        except EnvironmentError as e:
+            warnings.warn(str(e), UserWarning)
+        else:
+            for line in fo:
+                name_value = [s.strip() for s in line.split(':', 1)]
+                if len(name_value) != 2:
+                    continue
+                name, value = name_value
+                if not info or name in info[-1]:  # next processor
+                    info.append({})
+                info[-1][name] = value
+            fo.close()
+        self.__class__.info = info
+
+    def _not_impl(self):
+        pass
+
+    # Athlon
+
+    def _is_AMD(self):
+        return self.info[0]['vendor_id'] == 'AuthenticAMD'
+
+    def _is_AthlonK6_2(self):
+        return self._is_AMD() and self.info[0]['model'] == '2'
+
+    def _is_AthlonK6_3(self):
+        return self._is_AMD() and self.info[0]['model'] == '3'
+
+    def _is_AthlonK6(self):
+        return re.match(r'.*?AMD-K6', self.info[0]['model name']) is not None
+
+    def _is_AthlonK7(self):
+        return re.match(r'.*?AMD-K7', self.info[0]['model name']) is not None
+
+    def _is_AthlonMP(self):
+        return re.match(r'.*?Athlon\(tm\) MP\b',
+                        self.info[0]['model name']) is not None
+
+    def _is_AMD64(self):
+        return self.is_AMD() and self.info[0]['family'] == '15'
+
+    def _is_Athlon64(self):
+        return re.match(r'.*?Athlon\(tm\) 64\b',
+                        self.info[0]['model name']) is not None
+
+    def _is_AthlonHX(self):
+        return re.match(r'.*?Athlon HX\b',
+                        self.info[0]['model name']) is not None
+
+    def _is_Opteron(self):
+        return re.match(r'.*?Opteron\b',
+                        self.info[0]['model name']) is not None
+
+    def _is_Hammer(self):
+        return re.match(r'.*?Hammer\b',
+                        self.info[0]['model name']) is not None
+
+    # Alpha
+
+    def _is_Alpha(self):
+        return self.info[0]['cpu'] == 'Alpha'
+
+    def _is_EV4(self):
+        return self.is_Alpha() and self.info[0]['cpu model'] == 'EV4'
+
+    def _is_EV5(self):
+        return self.is_Alpha() and self.info[0]['cpu model'] == 'EV5'
+
+    def _is_EV56(self):
+        return self.is_Alpha() and self.info[0]['cpu model'] == 'EV56'
+
+    def _is_PCA56(self):
+        return self.is_Alpha() and self.info[0]['cpu model'] == 'PCA56'
+
+    # Intel
+
+    #XXX
+    _is_i386 = _not_impl
+
+    def _is_Intel(self):
+        return self.info[0]['vendor_id'] == 'GenuineIntel'
+
+    def _is_i486(self):
+        return self.info[0]['cpu'] == 'i486'
+
+    def _is_i586(self):
+        return self.is_Intel() and self.info[0]['cpu family'] == '5'
+
+    def _is_i686(self):
+        return self.is_Intel() and self.info[0]['cpu family'] == '6'
+
+    def _is_Celeron(self):
+        return re.match(r'.*?Celeron',
+                        self.info[0]['model name']) is not None
+
+    def _is_Pentium(self):
+        return re.match(r'.*?Pentium',
+                        self.info[0]['model name']) is not None
+
+    def _is_PentiumII(self):
+        return re.match(r'.*?Pentium.*?II\b',
+                        self.info[0]['model name']) is not None
+
+    def _is_PentiumPro(self):
+        return re.match(r'.*?PentiumPro\b',
+                        self.info[0]['model name']) is not None
+
+    def _is_PentiumMMX(self):
+        return re.match(r'.*?Pentium.*?MMX\b',
+                        self.info[0]['model name']) is not None
+
+    def _is_PentiumIII(self):
+        return re.match(r'.*?Pentium.*?III\b',
+                        self.info[0]['model name']) is not None
+
+    def _is_PentiumIV(self):
+        return re.match(r'.*?Pentium.*?(IV|4)\b',
+                        self.info[0]['model name']) is not None
+
+    def _is_PentiumM(self):
+        return re.match(r'.*?Pentium.*?M\b',
+                        self.info[0]['model name']) is not None
+
+    def _is_Prescott(self):
+        return self.is_PentiumIV() and self.has_sse3()
+
+    def _is_Nocona(self):
+        return (self.is_Intel() and
+                self.info[0]['cpu family'] in ('6', '15') and
+                # two s sse3; three s ssse3 not the same thing, this is fine
+                (self.has_sse3() and not self.has_ssse3()) and
+                re.match(r'.*?\blm\b', self.info[0]['flags']) is not None)
+
+    def _is_Core2(self):
+        return (self.is_64bit() and self.is_Intel() and
+                re.match(r'.*?Core\(TM\)2\b',
+                         self.info[0]['model name']) is not None)
+
+    def _is_Itanium(self):
+        return re.match(r'.*?Itanium\b',
+                        self.info[0]['family']) is not None
+
+    def _is_XEON(self):
+        return re.match(r'.*?XEON\b',
+                        self.info[0]['model name'], re.IGNORECASE) is not None
+
+    _is_Xeon = _is_XEON
+
+    # Power
+    def _is_Power(self):
+        return re.match(r'.*POWER.*',
+                       self.info[0]['cpu']) is not None
+
+    def _is_Power7(self):
+        return re.match(r'.*POWER7.*',
+                       self.info[0]['cpu']) is not None
+
+    def _is_Power8(self):
+        return re.match(r'.*POWER8.*',
+                       self.info[0]['cpu']) is not None
+
+    def _is_Power9(self):
+        return re.match(r'.*POWER9.*',
+                       self.info[0]['cpu']) is not None
+
+    def _has_Altivec(self):
+        return re.match(r'.*altivec\ supported.*',
+                       self.info[0]['cpu']) is not None
+
+    # Varia
+
+    def _is_singleCPU(self):
+        return len(self.info) == 1
+
+    def _getNCPUs(self):
+        return len(self.info)
+
+    def _has_fdiv_bug(self):
+        return self.info[0]['fdiv_bug'] == 'yes'
+
+    def _has_f00f_bug(self):
+        return self.info[0]['f00f_bug'] == 'yes'
+
+    def _has_mmx(self):
+        return re.match(r'.*?\bmmx\b', self.info[0]['flags']) is not None
+
+    def _has_sse(self):
+        return re.match(r'.*?\bsse\b', self.info[0]['flags']) is not None
+
+    def _has_sse2(self):
+        return re.match(r'.*?\bsse2\b', self.info[0]['flags']) is not None
+
+    def _has_sse3(self):
+        return re.match(r'.*?\bpni\b', self.info[0]['flags']) is not None
+
+    def _has_ssse3(self):
+        return re.match(r'.*?\bssse3\b', self.info[0]['flags']) is not None
+
+    def _has_3dnow(self):
+        return re.match(r'.*?\b3dnow\b', self.info[0]['flags']) is not None
+
+    def _has_3dnowext(self):
+        return re.match(r'.*?\b3dnowext\b', self.info[0]['flags']) is not None
+
+
+class IRIXCPUInfo(CPUInfoBase):
+    info = None
+
+    def __init__(self):
+        if self.info is not None:
+            return
+        info = key_value_from_command('sysconf', sep=' ',
+                                      successful_status=(0, 1))
+        self.__class__.info = info
+
+    def _not_impl(self):
+        pass
+
+    def _is_singleCPU(self):
+        return self.info.get('NUM_PROCESSORS') == '1'
+
+    def _getNCPUs(self):
+        return int(self.info.get('NUM_PROCESSORS', 1))
+
+    def __cputype(self, n):
+        return self.info.get('PROCESSORS').split()[0].lower() == 'r%s' % (n)
+
+    def _is_r2000(self):
+        return self.__cputype(2000)
+
+    def _is_r3000(self):
+        return self.__cputype(3000)
+
+    def _is_r3900(self):
+        return self.__cputype(3900)
+
+    def _is_r4000(self):
+        return self.__cputype(4000)
+
+    def _is_r4100(self):
+        return self.__cputype(4100)
+
+    def _is_r4300(self):
+        return self.__cputype(4300)
+
+    def _is_r4400(self):
+        return self.__cputype(4400)
+
+    def _is_r4600(self):
+        return self.__cputype(4600)
+
+    def _is_r4650(self):
+        return self.__cputype(4650)
+
+    def _is_r5000(self):
+        return self.__cputype(5000)
+
+    def _is_r6000(self):
+        return self.__cputype(6000)
+
+    def _is_r8000(self):
+        return self.__cputype(8000)
+
+    def _is_r10000(self):
+        return self.__cputype(10000)
+
+    def _is_r12000(self):
+        return self.__cputype(12000)
+
+    def _is_rorion(self):
+        return self.__cputype('orion')
+
+    def get_ip(self):
+        try:
+            return self.info.get('MACHINE')
+        except:
+            pass
+
+    def __machine(self, n):
+        return self.info.get('MACHINE').lower() == 'ip%s' % (n)
+
+    def _is_IP19(self):
+        return self.__machine(19)
+
+    def _is_IP20(self):
+        return self.__machine(20)
+
+    def _is_IP21(self):
+        return self.__machine(21)
+
+    def _is_IP22(self):
+        return self.__machine(22)
+
+    def _is_IP22_4k(self):
+        return self.__machine(22) and self._is_r4000()
+
+    def _is_IP22_5k(self):
+        return self.__machine(22) and self._is_r5000()
+
+    def _is_IP24(self):
+        return self.__machine(24)
+
+    def _is_IP25(self):
+        return self.__machine(25)
+
+    def _is_IP26(self):
+        return self.__machine(26)
+
+    def _is_IP27(self):
+        return self.__machine(27)
+
+    def _is_IP28(self):
+        return self.__machine(28)
+
+    def _is_IP30(self):
+        return self.__machine(30)
+
+    def _is_IP32(self):
+        return self.__machine(32)
+
+    def _is_IP32_5k(self):
+        return self.__machine(32) and self._is_r5000()
+
+    def _is_IP32_10k(self):
+        return self.__machine(32) and self._is_r10000()
+
+
+class DarwinCPUInfo(CPUInfoBase):
+    info = None
+
+    def __init__(self):
+        if self.info is not None:
+            return
+        info = command_info(arch='arch',
+                            machine='machine')
+        info['sysctl_hw'] = key_value_from_command(['sysctl', 'hw'], sep='=')
+        self.__class__.info = info
+
+    def _not_impl(self): pass
+
+    def _getNCPUs(self):
+        return int(self.info['sysctl_hw'].get('hw.ncpu', 1))
+
+    def _is_Power_Macintosh(self):
+        return self.info['sysctl_hw']['hw.machine'] == 'Power Macintosh'
+
+    def _is_i386(self):
+        return self.info['arch'] == 'i386'
+
+    def _is_ppc(self):
+        return self.info['arch'] == 'ppc'
+
+    def __machine(self, n):
+        return self.info['machine'] == 'ppc%s' % n
+
+    def _is_ppc601(self): return self.__machine(601)
+
+    def _is_ppc602(self): return self.__machine(602)
+
+    def _is_ppc603(self): return self.__machine(603)
+
+    def _is_ppc603e(self): return self.__machine('603e')
+
+    def _is_ppc604(self): return self.__machine(604)
+
+    def _is_ppc604e(self): return self.__machine('604e')
+
+    def _is_ppc620(self): return self.__machine(620)
+
+    def _is_ppc630(self): return self.__machine(630)
+
+    def _is_ppc740(self): return self.__machine(740)
+
+    def _is_ppc7400(self): return self.__machine(7400)
+
+    def _is_ppc7450(self): return self.__machine(7450)
+
+    def _is_ppc750(self): return self.__machine(750)
+
+    def _is_ppc403(self): return self.__machine(403)
+
+    def _is_ppc505(self): return self.__machine(505)
+
+    def _is_ppc801(self): return self.__machine(801)
+
+    def _is_ppc821(self): return self.__machine(821)
+
+    def _is_ppc823(self): return self.__machine(823)
+
+    def _is_ppc860(self): return self.__machine(860)
+
+class NetBSDCPUInfo(CPUInfoBase):
+	info = None
+
+	def __init__(self):
+		if self.info is not None:
+			return
+		info = {}
+		info['sysctl_hw'] = key_value_from_command(['sysctl', 'hw'], sep='=')
+		info['arch'] = info['sysctl_hw'].get('hw.machine_arch', 1)
+		info['machine'] = info['sysctl_hw'].get('hw.machine', 1)
+		self.__class__.info = info
+
+	def _not_impl(self): pass
+
+	def _getNCPUs(self):
+		return int(self.info['sysctl_hw'].get('hw.ncpu', 1))
+
+	def _is_Intel(self):
+		if self.info['sysctl_hw'].get('hw.model', "")[0:5] == 'Intel':
+			return True
+		return False
+
+	def _is_AMD(self):
+		if self.info['sysctl_hw'].get('hw.model', "")[0:3] == 'AMD':
+			return True
+		return False
+
+class SunOSCPUInfo(CPUInfoBase):
+    info = None
+
+    def __init__(self):
+        if self.info is not None:
+            return
+        info = command_info(arch='arch',
+                            mach='mach',
+                            uname_i=['uname', '-i'],
+                            isainfo_b=['isainfo', '-b'],
+                            isainfo_n=['isainfo', '-n'],
+        )
+        info['uname_X'] = key_value_from_command(['uname', '-X'], sep='=')
+        for line in command_by_line(['psrinfo', '-v', '0']):
+            m = re.match(r'\s*The (?P<p>[\w\d]+) processor operates at', line)
+            if m:
+                info['processor'] = m.group('p')
+                break
+        self.__class__.info = info
+
+    def _not_impl(self):
+        pass
+
+    def _is_i386(self):
+        return self.info['isainfo_n'] == 'i386'
+
+    def _is_sparc(self):
+        return self.info['isainfo_n'] == 'sparc'
+
+    def _is_sparcv9(self):
+        return self.info['isainfo_n'] == 'sparcv9'
+
+    def _getNCPUs(self):
+        return int(self.info['uname_X'].get('NumCPU', 1))
+
+    def _is_sun4(self):
+        return self.info['arch'] == 'sun4'
+
+    def _is_SUNW(self):
+        return re.match(r'SUNW', self.info['uname_i']) is not None
+
+    def _is_sparcstation5(self):
+        return re.match(r'.*SPARCstation-5', self.info['uname_i']) is not None
+
+    def _is_ultra1(self):
+        return re.match(r'.*Ultra-1', self.info['uname_i']) is not None
+
+    def _is_ultra250(self):
+        return re.match(r'.*Ultra-250', self.info['uname_i']) is not None
+
+    def _is_ultra2(self):
+        return re.match(r'.*Ultra-2', self.info['uname_i']) is not None
+
+    def _is_ultra30(self):
+        return re.match(r'.*Ultra-30', self.info['uname_i']) is not None
+
+    def _is_ultra4(self):
+        return re.match(r'.*Ultra-4', self.info['uname_i']) is not None
+
+    def _is_ultra5_10(self):
+        return re.match(r'.*Ultra-5_10', self.info['uname_i']) is not None
+
+    def _is_ultra5(self):
+        return re.match(r'.*Ultra-5', self.info['uname_i']) is not None
+
+    def _is_ultra60(self):
+        return re.match(r'.*Ultra-60', self.info['uname_i']) is not None
+
+    def _is_ultra80(self):
+        return re.match(r'.*Ultra-80', self.info['uname_i']) is not None
+
+    def _is_ultraenterprice(self):
+        return re.match(r'.*Ultra-Enterprise', self.info['uname_i']) is not None
+
+    def _is_ultraenterprice10k(self):
+        return re.match(r'.*Ultra-Enterprise-10000', self.info['uname_i']) is not None
+
+    def _is_sunfire(self):
+        return re.match(r'.*Sun-Fire', self.info['uname_i']) is not None
+
+    def _is_ultra(self):
+        return re.match(r'.*Ultra', self.info['uname_i']) is not None
+
+    def _is_cpusparcv7(self):
+        return self.info['processor'] == 'sparcv7'
+
+    def _is_cpusparcv8(self):
+        return self.info['processor'] == 'sparcv8'
+
+    def _is_cpusparcv9(self):
+        return self.info['processor'] == 'sparcv9'
+
+
+class Win32CPUInfo(CPUInfoBase):
+    info = None
+    pkey = r"HARDWARE\DESCRIPTION\System\CentralProcessor"
+    # XXX: what does the value of
+    #   HKEY_LOCAL_MACHINE\HARDWARE\DESCRIPTION\System\CentralProcessor\0
+    # mean?
+
+    def __init__(self):
+        try:
+            import _winreg
+        except ImportError:  # Python 3
+            import winreg as _winreg
+
+        if self.info is not None:
+            return
+        info = []
+        try:
+            #XXX: Bad style to use so long `try:...except:...`. Fix it!
+
+            prgx = re.compile(r"family\s+(?P<FML>\d+)\s+model\s+(?P<MDL>\d+)"
+                              r"\s+stepping\s+(?P<STP>\d+)", re.IGNORECASE)
+            chnd = _winreg.OpenKey(_winreg.HKEY_LOCAL_MACHINE, self.pkey)
+            pnum = 0
+            while 1:
+                try:
+                    proc = _winreg.EnumKey(chnd, pnum)
+                except _winreg.error:
+                    break
+                else:
+                    pnum += 1
+                    info.append({"Processor": proc})
+                    phnd = _winreg.OpenKey(chnd, proc)
+                    pidx = 0
+                    while True:
+                        try:
+                            name, value, vtpe = _winreg.EnumValue(phnd, pidx)
+                        except _winreg.error:
+                            break
+                        else:
+                            pidx = pidx + 1
+                            info[-1][name] = value
+                            if name == "Identifier":
+                                srch = prgx.search(value)
+                                if srch:
+                                    info[-1]["Family"] = int(srch.group("FML"))
+                                    info[-1]["Model"] = int(srch.group("MDL"))
+                                    info[-1]["Stepping"] = int(srch.group("STP"))
+        except:
+            print(sys.exc_value, '(ignoring)')
+        self.__class__.info = info
+
+    def _not_impl(self):
+        pass
+
+    # Athlon
+
+    def _is_AMD(self):
+        return self.info[0]['VendorIdentifier'] == 'AuthenticAMD'
+
+    def _is_Am486(self):
+        return self.is_AMD() and self.info[0]['Family'] == 4
+
+    def _is_Am5x86(self):
+        return self.is_AMD() and self.info[0]['Family'] == 4
+
+    def _is_AMDK5(self):
+        return (self.is_AMD() and self.info[0]['Family'] == 5 and
+                self.info[0]['Model'] in [0, 1, 2, 3])
+
+    def _is_AMDK6(self):
+        return (self.is_AMD() and self.info[0]['Family'] == 5 and
+                self.info[0]['Model'] in [6, 7])
+
+    def _is_AMDK6_2(self):
+        return (self.is_AMD() and self.info[0]['Family'] == 5 and
+                self.info[0]['Model'] == 8)
+
+    def _is_AMDK6_3(self):
+        return (self.is_AMD() and self.info[0]['Family'] == 5 and
+                self.info[0]['Model'] == 9)
+
+    def _is_AMDK7(self):
+        return self.is_AMD() and self.info[0]['Family'] == 6
+
+    # To reliably distinguish between the different types of AMD64 chips
+    # (Athlon64, Operton, Athlon64 X2, Semperon, Turion 64, etc.) would
+    # require looking at the 'brand' from cpuid
+
+    def _is_AMD64(self):
+        return self.is_AMD() and self.info[0]['Family'] == 15
+
+    # Intel
+
+    def _is_Intel(self):
+        return self.info[0]['VendorIdentifier'] == 'GenuineIntel'
+
+    def _is_i386(self):
+        return self.info[0]['Family'] == 3
+
+    def _is_i486(self):
+        return self.info[0]['Family'] == 4
+
+    def _is_i586(self):
+        return self.is_Intel() and self.info[0]['Family'] == 5
+
+    def _is_i686(self):
+        return self.is_Intel() and self.info[0]['Family'] == 6
+
+    def _is_Pentium(self):
+        return self.is_Intel() and self.info[0]['Family'] == 5
+
+    def _is_PentiumMMX(self):
+        return (self.is_Intel() and self.info[0]['Family'] == 5 and
+                self.info[0]['Model'] == 4)
+
+    def _is_PentiumPro(self):
+        return (self.is_Intel() and self.info[0]['Family'] == 6 and
+                self.info[0]['Model'] == 1)
+
+    def _is_PentiumII(self):
+        return (self.is_Intel() and self.info[0]['Family'] == 6 and
+                self.info[0]['Model'] in [3, 5, 6])
+
+    def _is_PentiumIII(self):
+        return (self.is_Intel() and self.info[0]['Family'] == 6 and
+                self.info[0]['Model'] in [7, 8, 9, 10, 11])
+
+    def _is_PentiumIV(self):
+        return self.is_Intel() and self.info[0]['Family'] == 15
+
+    def _is_PentiumM(self):
+        return (self.is_Intel() and self.info[0]['Family'] == 6 and
+                self.info[0]['Model'] in [9, 13, 14])
+
+    def _is_Core2(self):
+        return (self.is_Intel() and self.info[0]['Family'] == 6 and
+                self.info[0]['Model'] in [15, 16, 17])
+
+    # Varia
+
+    def _is_singleCPU(self):
+        return len(self.info) == 1
+
+    def _getNCPUs(self):
+        return len(self.info)
+
+    def _has_mmx(self):
+        if self.is_Intel():
+            return ((self.info[0]['Family'] == 5 and
+                     self.info[0]['Model'] == 4) or
+                    (self.info[0]['Family'] in [6, 15]))
+        elif self.is_AMD():
+            return self.info[0]['Family'] in [5, 6, 15]
+        else:
+            return False
+
+    def _has_sse(self):
+        if self.is_Intel():
+            return ((self.info[0]['Family'] == 6 and
+                     self.info[0]['Model'] in [7, 8, 9, 10, 11]) or
+                    self.info[0]['Family'] == 15)
+        elif self.is_AMD():
+            return ((self.info[0]['Family'] == 6 and
+                     self.info[0]['Model'] in [6, 7, 8, 10]) or
+                    self.info[0]['Family'] == 15)
+        else:
+            return False
+
+    def _has_sse2(self):
+        if self.is_Intel():
+            return self.is_Pentium4() or self.is_PentiumM() or self.is_Core2()
+        elif self.is_AMD():
+            return self.is_AMD64()
+        else:
+            return False
+
+    def _has_3dnow(self):
+        return self.is_AMD() and self.info[0]['Family'] in [5, 6, 15]
+
+    def _has_3dnowext(self):
+        return self.is_AMD() and self.info[0]['Family'] in [6, 15]
+
+
+if sys.platform.startswith('linux'):  # variations: linux2,linux-i386 (any others?)
+    cpuinfo = LinuxCPUInfo
+elif sys.platform.startswith('irix'):
+    cpuinfo = IRIXCPUInfo
+elif sys.platform == 'darwin':
+    cpuinfo = DarwinCPUInfo
+elif sys.platform[0:6] == 'netbsd':
+    cpuinfo = NetBSDCPUInfo
+elif sys.platform.startswith('sunos'):
+    cpuinfo = SunOSCPUInfo
+elif sys.platform.startswith('win32'):
+    cpuinfo = Win32CPUInfo
+elif sys.platform.startswith('cygwin'):
+    cpuinfo = LinuxCPUInfo
+#XXX: other OS's. Eg. use _winreg on Win32. Or os.uname on unices.
+else:
+    cpuinfo = CPUInfoBase
+
+cpu = cpuinfo()
+
+if __name__ == "__main__":
+
+    cpu.is_blaa()
+    cpu.is_Intel()
+    cpu.is_Alpha()
+
+    info = []
+    for name in dir(cpuinfo):
+        if name[0] == '_' and name[1] != '_':
+            r = getattr(cpu, name[1:])()
+            if r:
+                if r != 1:
+                    info.append('%s=%s' % (name[1:], r))
+                else:
+                    info.append(name[1:])
+    print('CPU information: ' + ' '.join(info))

numexpr/expressions.py

@@ -0,0 +1,546 @@
+###################################################################
+#  Numexpr - Fast numerical array expression evaluator for NumPy.
+#
+#      License: MIT
+#      Author:  See AUTHORS.txt
+#
+#  See LICENSE.txt and LICENSES/*.txt for details about copyright and
+#  rights to use.
+####################################################################
+
+__all__ = ['E']
+
+import operator
+import sys
+import threading
+
+import numpy
+
+# Declare a double type that does not exist in Python space
+double = numpy.double
+
+# The default kind for undeclared variables
+default_kind = 'double'
+int_ = numpy.int32
+long_ = numpy.int64
+
+type_to_kind = {bool: 'bool', int_: 'int', long_: 'long', float: 'float',
+                double: 'double', complex: 'complex', bytes: 'bytes', str: 'str'}
+kind_to_type = {'bool': bool, 'int': int_, 'long': long_, 'float': float,
+                'double': double, 'complex': complex, 'bytes': bytes, 'str': str}
+kind_rank = ('bool', 'int', 'long', 'float', 'double', 'complex', 'none')
+scalar_constant_types = [bool, int_, int, float, double, complex, bytes, str]
+
+scalar_constant_types = tuple(scalar_constant_types)
+
+from numexpr import interpreter
+
+
+class Expression():
+
+    def __getattr__(self, name):
+        if name.startswith('_'):
+            try:
+                return self.__dict__[name]
+            except KeyError:
+                raise AttributeError
+        else:
+            return VariableNode(name, default_kind)
+
+
+E = Expression()
+
+
+class Context(threading.local):
+
+    def get(self, value, default):
+        return self.__dict__.get(value, default)
+
+    def get_current_context(self):
+        return self.__dict__
+
+    def set_new_context(self, dict_):
+        self.__dict__.update(dict_)
+
+# This will be called each time the local object is used in a separate thread
+_context = Context()
+
+
+def get_optimization():
+    return _context.get('optimization', 'none')
+
+
+# helper functions for creating __magic__ methods
+def ophelper(f):
+    def func(*args):
+        args = list(args)
+        for i, x in enumerate(args):
+            if isConstant(x):
+                args[i] = x = ConstantNode(x)
+            if not isinstance(x, ExpressionNode):
+                raise TypeError("unsupported object type: %s" % type(x))
+        return f(*args)
+
+    func.__name__ = f.__name__
+    func.__doc__ = f.__doc__
+    func.__dict__.update(f.__dict__)
+    return func
+
+
+def allConstantNodes(args):
+    "returns True if args are all ConstantNodes."
+    for x in args:
+        if not isinstance(x, ConstantNode):
+            return False
+    return True
+
+
+def isConstant(ex):
+    "Returns True if ex is a constant scalar of an allowed type."
+    return isinstance(ex, scalar_constant_types)
+
+
+def commonKind(nodes):
+    node_kinds = [node.astKind for node in nodes]
+    str_count = node_kinds.count('bytes') + node_kinds.count('str')
+    if 0 < str_count < len(node_kinds):  # some args are strings, but not all
+        raise TypeError("strings can only be operated with strings")
+    if str_count > 0:  # if there are some, all of them must be
+        return 'bytes'
+    n = -1
+    for x in nodes:
+        n = max(n, kind_rank.index(x.astKind))
+    return kind_rank[n]
+
+
+max_int32 = 2147483647
+min_int32 = -max_int32 - 1
+
+
+def bestConstantType(x):
+    # ``numpy.string_`` is a subclass of ``bytes``
+    if isinstance(x, (bytes, str)):
+        return bytes
+    # Numeric conversion to boolean values is not tried because
+    # ``bool(1) == True`` (same for 0 and False), so 0 and 1 would be
+    # interpreted as booleans when ``False`` and ``True`` are already
+    # supported.
+    if isinstance(x, (bool, numpy.bool_)):
+        return bool
+    # ``long`` objects are kept as is to allow the user to force
+    # promotion of results by using long constants, e.g. by operating
+    # a 32-bit array with a long (64-bit) constant.
+    if isinstance(x, (long_, numpy.int64)):
+        return long_
+    # ``double`` objects are kept as is to allow the user to force
+    # promotion of results by using double constants, e.g. by operating
+    # a float (32-bit) array with a double (64-bit) constant.
+    if isinstance(x, double):
+        return double
+    if isinstance(x, numpy.float32):
+        return float
+    if isinstance(x, (int, numpy.integer)):
+        # Constants needing more than 32 bits are always
+        # considered ``long``, *regardless of the platform*, so we
+        # can clearly tell 32- and 64-bit constants apart.
+        if not (min_int32 <= x <= max_int32):
+            return long_
+        return int_
+    # The duality of float and double in Python avoids that we have to list
+    # ``double`` too.
+    for converter in float, complex:
+        try:
+            y = converter(x)
+        except Exception as err:
+            continue
+        if y == x or numpy.isnan(y):
+            return converter
+
+
+def getKind(x):
+    converter = bestConstantType(x)
+    return type_to_kind[converter]
+
+
+def binop(opname, reversed=False, kind=None):
+    # Getting the named method from self (after reversal) does not
+    # always work (e.g. int constants do not have a __lt__ method).
+    opfunc = getattr(operator, "__%s__" % opname)
+
+    @ophelper
+    def operation(self, other):
+        if reversed:
+            self, other = other, self
+        if allConstantNodes([self, other]):
+            return ConstantNode(opfunc(self.value, other.value))
+        else:
+            return OpNode(opname, (self, other), kind=kind)
+
+    return operation
+
+
+def func(func, minkind=None, maxkind=None):
+    @ophelper
+    def function(*args):
+        if allConstantNodes(args):
+            return ConstantNode(func(*[x.value for x in args]))
+        kind = commonKind(args)
+        if kind in ('int', 'long'):
+            if func.__name__ not in ('copy', 'abs', 'ones_like', 'round', 'sign'):
+                # except for these special functions (which return ints for int inputs in NumPy)
+                # just do a cast to double
+                # FIXME: 'fmod' outputs ints for NumPy when inputs are ints, but need to
+                # add new function signatures FUNC_LLL FUNC_III to support this
+                kind = 'double'
+        else:
+            # Apply regular casting rules
+            if minkind and kind_rank.index(minkind) > kind_rank.index(kind):
+                kind = minkind
+            if maxkind and kind_rank.index(maxkind) < kind_rank.index(kind):
+                kind = maxkind
+        return FuncNode(func.__name__, args, kind)
+
+    return function
+
+
+@ophelper
+def where_func(a, b, c):
+    if isinstance(a, ConstantNode):
+        return b if a.value else c
+    if allConstantNodes([a, b, c]):
+        return ConstantNode(numpy.where(a, b, c))
+    return FuncNode('where', [a, b, c])
+
+
+def encode_axis(axis):
+    if isinstance(axis, ConstantNode):
+        axis = axis.value
+    if axis is None:
+        axis = interpreter.allaxes
+    else:
+        if axis < 0:
+            raise ValueError("negative axis are not supported")
+        if axis > 254:
+            raise ValueError("cannot encode axis")
+    return RawNode(axis)
+
+
+def gen_reduce_axis_func(name):
+    def _func(a, axis=None):
+        axis = encode_axis(axis)
+        if isinstance(a, ConstantNode):
+            return a
+        if isinstance(a, (bool, int_, long_, float, double, complex)):
+            a = ConstantNode(a)
+        return FuncNode(name, [a, axis], kind=a.astKind)
+    return _func
+
+
+@ophelper
+def contains_func(a, b):
+    return FuncNode('contains', [a, b], kind='bool')
+
+
+@ophelper
+def div_op(a, b):
+    if get_optimization() in ('moderate', 'aggressive'):
+        if (isinstance(b, ConstantNode) and
+                (a.astKind == b.astKind) and
+                    a.astKind in ('float', 'double', 'complex')):
+            return OpNode('mul', [a, ConstantNode(1. / b.value)])
+    return OpNode('div', [a, b])
+
+
+@ophelper
+def truediv_op(a, b):
+    if get_optimization() in ('moderate', 'aggressive'):
+        if (isinstance(b, ConstantNode) and
+                (a.astKind == b.astKind) and
+                    a.astKind in ('float', 'double', 'complex')):
+            return OpNode('mul', [a, ConstantNode(1. / b.value)])
+    kind = commonKind([a, b])
+    if kind in ('bool', 'int', 'long'):
+        kind = 'double'
+    return OpNode('div', [a, b], kind=kind)
+
+
+@ophelper
+def rtruediv_op(a, b):
+    return truediv_op(b, a)
+
+
+@ophelper
+def pow_op(a, b):
+
+    if isinstance(b, ConstantNode):
+        x = b.value
+        if (    a.astKind in ('int', 'long') and
+                b.astKind in ('int', 'long') and x < 0) :
+            raise ValueError(
+                'Integers to negative integer powers are not allowed.')
+        if get_optimization() == 'aggressive':
+            RANGE = 50  # Approximate break even point with pow(x,y)
+            # Optimize all integral and half integral powers in [-RANGE, RANGE]
+            # Note: for complex numbers RANGE could be larger.
+            if (int(2 * x) == 2 * x) and (-RANGE <= abs(x) <= RANGE):
+                n = int_(abs(x))
+                ishalfpower = int_(abs(2 * x)) % 2
+
+                def multiply(x, y):
+                    if x is None: return y
+                    return OpNode('mul', [x, y])
+
+                r = None
+                p = a
+                mask = 1
+                while True:
+                    if (n & mask):
+                        r = multiply(r, p)
+                    mask <<= 1
+                    if mask > n:
+                        break
+                    p = OpNode('mul', [p, p])
+                if ishalfpower:
+                    kind = commonKind([a])
+                    if kind in ('int', 'long'):
+                        kind = 'double'
+                    r = multiply(r, OpNode('sqrt', [a], kind))
+                if r is None:
+                    r = OpNode('ones_like', [a])
+                if x < 0:
+                    # Issue #428
+                    r = truediv_op(ConstantNode(1), r)
+                return r
+        if get_optimization() in ('moderate', 'aggressive'):
+            if x == -1:
+                return OpNode('div', [ConstantNode(1), a])
+            if x == 0:
+                return OpNode('ones_like', [a])
+            if x == 0.5:
+                kind = a.astKind
+                if kind in ('int', 'long'): kind = 'double'
+                return FuncNode('sqrt', [a], kind=kind)
+            if x == 1:
+                return a
+            if x == 2:
+                return OpNode('mul', [a, a])
+    return OpNode('pow', [a, b])
+
+# The functions and the minimum and maximum types accepted
+numpy.expm1x = numpy.expm1
+functions = {
+    'copy': func(numpy.copy),
+    'ones_like': func(numpy.ones_like),
+    'sqrt': func(numpy.sqrt, 'float'),
+
+    'sin': func(numpy.sin, 'float'),
+    'cos': func(numpy.cos, 'float'),
+    'tan': func(numpy.tan, 'float'),
+    'arcsin': func(numpy.arcsin, 'float'),
+    'arccos': func(numpy.arccos, 'float'),
+    'arctan': func(numpy.arctan, 'float'),
+
+    'sinh': func(numpy.sinh, 'float'),
+    'cosh': func(numpy.cosh, 'float'),
+    'tanh': func(numpy.tanh, 'float'),
+    'arcsinh': func(numpy.arcsinh, 'float'),
+    'arccosh': func(numpy.arccosh, 'float'),
+    'arctanh': func(numpy.arctanh, 'float'),
+
+    'fmod': func(numpy.fmod, 'float'),
+    'arctan2': func(numpy.arctan2, 'float'),
+    'hypot': func(numpy.hypot, 'double'),
+    'nextafter': func(numpy.nextafter, 'double'),
+    'copysign': func(numpy.copysign, 'double'),
+    'maximum': func(numpy.maximum, 'double'),
+    'minimum': func(numpy.minimum, 'double'),
+
+
+    'log': func(numpy.log, 'float'),
+    'log1p': func(numpy.log1p, 'float'),
+    'log10': func(numpy.log10, 'float'),
+    'log2': func(numpy.log2, 'float'),
+    'exp': func(numpy.exp, 'float'),
+    'expm1': func(numpy.expm1, 'float'),
+
+    'abs': func(numpy.absolute, 'float'),
+    'ceil': func(numpy.ceil, 'float', 'double'),
+    'floor': func(numpy.floor, 'float', 'double'),
+    'round': func(numpy.round, 'double'),
+    'trunc': func(numpy.trunc, 'double'),
+    'sign': func(numpy.sign, 'double'),
+
+    'where': where_func,
+
+    'real': func(numpy.real, 'double', 'double'),
+    'imag': func(numpy.imag, 'double', 'double'),
+    'complex': func(complex, 'complex'),
+    'conj': func(numpy.conj, 'complex'),
+
+    'isnan': func(numpy.isnan, 'double'),
+    'isfinite': func(numpy.isfinite, 'double'),
+    'isinf': func(numpy.isinf, 'double'),
+    'signbit': func(numpy.signbit, 'double'),
+
+    'sum': gen_reduce_axis_func('sum'),
+    'prod': gen_reduce_axis_func('prod'),
+    'min': gen_reduce_axis_func('min'),
+    'max': gen_reduce_axis_func('max'),
+    'contains': contains_func,
+}
+
+
+class ExpressionNode():
+    """
+    An object that represents a generic number object.
+
+    This implements the number special methods so that we can keep
+    track of how this object has been used.
+    """
+    astType = 'generic'
+
+    def __init__(self, value=None, kind=None, children=None):
+        self.value = value
+        if kind is None:
+            kind = 'none'
+        self.astKind = kind
+        if children is None:
+            self.children = ()
+        else:
+            self.children = tuple(children)
+
+    def get_real(self):
+        if self.astType == 'constant':
+            return ConstantNode(complex(self.value).real)
+        return OpNode('real', (self,), 'double')
+
+    real = property(get_real)
+
+    def get_imag(self):
+        if self.astType == 'constant':
+            return ConstantNode(complex(self.value).imag)
+        return OpNode('imag', (self,), 'double')
+
+    imag = property(get_imag)
+
+    def __str__(self):
+        return '%s(%s, %s, %s)' % (self.__class__.__name__, self.value,
+                                   self.astKind, self.children)
+
+    def __repr__(self):
+        return self.__str__()
+
+    def __neg__(self):
+        return OpNode('neg', (self,))
+
+    def __invert__(self):
+        return OpNode('invert', (self,))
+
+    def __pos__(self):
+        return self
+
+    # The next check is commented out. See #24 for more info.
+
+    def __bool__(self):
+        raise TypeError("You can't use Python's standard boolean operators in "
+                        "NumExpr expressions. You should use their bitwise "
+                        "counterparts instead: '&' instead of 'and', "
+                        "'|' instead of 'or', and '~' instead of 'not'.")
+
+    __add__ = __radd__ = binop('add')
+    __sub__ = binop('sub')
+    __rsub__ = binop('sub', reversed=True)
+    __mul__ = __rmul__ = binop('mul')
+    __truediv__ = truediv_op
+    __rtruediv__ = rtruediv_op
+    __floordiv__ = binop("floordiv")
+    __pow__ = pow_op
+    __rpow__ = binop('pow', reversed=True)
+    __mod__ = binop('mod')
+    __rmod__ = binop('mod', reversed=True)
+
+    __lshift__ = binop('lshift')
+    __rlshift__ = binop('lshift', reversed=True)
+    __rshift__ = binop('rshift')
+    __rrshift__ = binop('rshift', reversed=True)
+
+    # bitwise or logical operations
+    __and__ = binop('and')
+    __or__ = binop('or')
+    __xor__ = binop('xor')
+
+    __gt__ = binop('gt', kind='bool')
+    __ge__ = binop('ge', kind='bool')
+    __eq__ = binop('eq', kind='bool')
+    __ne__ = binop('ne', kind='bool')
+    __lt__ = binop('gt', reversed=True, kind='bool')
+    __le__ = binop('ge', reversed=True, kind='bool')
+
+
+class LeafNode(ExpressionNode):
+    leafNode = True
+
+
+class VariableNode(LeafNode):
+    astType = 'variable'
+
+    def __init__(self, value=None, kind=None, children=None):
+        LeafNode.__init__(self, value=value, kind=kind)
+
+
+class RawNode():
+    """
+    Used to pass raw integers to interpreter.
+    For instance, for selecting what function to use in func1.
+    Purposely don't inherit from ExpressionNode, since we don't wan't
+    this to be used for anything but being walked.
+    """
+    astType = 'raw'
+    astKind = 'none'
+
+    def __init__(self, value):
+        self.value = value
+        self.children = ()
+
+    def __str__(self):
+        return 'RawNode(%s)' % (self.value,)
+
+    __repr__ = __str__
+
+
+class ConstantNode(LeafNode):
+    astType = 'constant'
+
+    def __init__(self, value=None, children=None):
+        kind = getKind(value)
+        # Python float constants are double precision by default
+        if kind == 'float' and isinstance(value, float):
+            kind = 'double'
+        LeafNode.__init__(self, value=value, kind=kind)
+
+    def __neg__(self):
+        return ConstantNode(-self.value)
+
+    def __invert__(self):
+        return ConstantNode(~self.value)
+
+
+class OpNode(ExpressionNode):
+    astType = 'op'
+
+    def __init__(self, opcode=None, args=None, kind=None):
+        if (kind is None) and (args is not None):
+            kind = commonKind(args)
+        if kind=='bool': # handle bool*bool and bool+bool cases
+            opcode = 'and' if opcode=='mul' else opcode
+            opcode = 'or' if opcode=='add' else opcode
+        ExpressionNode.__init__(self, value=opcode, kind=kind, children=args)
+
+
+class FuncNode(OpNode):
+    def __init__(self, opcode=None, args=None, kind=None):
+        if (kind is None) and (args is not None):
+            kind = commonKind(args)
+        if opcode in ("isnan", "isfinite", "isinf", "signbit"): # bodge for boolean return functions
+            kind = 'bool'
+        OpNode.__init__(self, opcode, args, kind)

numexpr/functions.hpp

@@ -0,0 +1,235 @@
+// -*- c-mode -*-
+/*********************************************************************
+  Numexpr - Fast numerical array expression evaluator for NumPy.
+
+      License: MIT
+      Author:  See AUTHORS.txt
+
+  See LICENSE.txt for details about copyright and rights to use.
+**********************************************************************/
+
+/* These #if blocks make it easier to query this file, without having
+   to define every row function before #including it. */
+#ifndef FUNC_FF
+#define ELIDE_FUNC_FF
+#define FUNC_FF(...)
+#endif
+FUNC_FF(FUNC_SQRT_FF,    "sqrt_ff",     sqrtf,  sqrtf2,  vsSqrt)
+FUNC_FF(FUNC_SIN_FF,     "sin_ff",      sinf,   sinf2,   vsSin)
+FUNC_FF(FUNC_COS_FF,     "cos_ff",      cosf,   cosf2,   vsCos)
+FUNC_FF(FUNC_TAN_FF,     "tan_ff",      tanf,   tanf2,   vsTan)
+FUNC_FF(FUNC_ARCSIN_FF,  "arcsin_ff",   asinf,  asinf2,  vsAsin)
+FUNC_FF(FUNC_ARCCOS_FF,  "arccos_ff",   acosf,  acosf2,  vsAcos)
+FUNC_FF(FUNC_ARCTAN_FF,  "arctan_ff",   atanf,  atanf2,  vsAtan)
+FUNC_FF(FUNC_SINH_FF,    "sinh_ff",     sinhf,  sinhf2,  vsSinh)
+FUNC_FF(FUNC_COSH_FF,    "cosh_ff",     coshf,  coshf2,  vsCosh)
+FUNC_FF(FUNC_TANH_FF,    "tanh_ff",     tanhf,  tanhf2,  vsTanh)
+FUNC_FF(FUNC_ARCSINH_FF, "arcsinh_ff",  asinhf, asinhf2, vsAsinh)
+FUNC_FF(FUNC_ARCCOSH_FF, "arccosh_ff",  acoshf, acoshf2, vsAcosh)
+FUNC_FF(FUNC_ARCTANH_FF, "arctanh_ff",  atanhf, atanhf2, vsAtanh)
+FUNC_FF(FUNC_LOG_FF,     "log_ff",      logf,   logf2,   vsLn)
+FUNC_FF(FUNC_LOG1P_FF,   "log1p_ff",    log1pf, log1pf2, vsLog1p)
+FUNC_FF(FUNC_LOG10_FF,   "log10_ff",    log10f, log10f2, vsLog10)
+FUNC_FF(FUNC_LOG2_FF,    "log2_ff",     log2f,  log2f2,  vsLog2)
+FUNC_FF(FUNC_EXP_FF,     "exp_ff",      expf,   expf2,   vsExp)
+FUNC_FF(FUNC_EXPM1_FF,   "expm1_ff",    expm1f, expm1f2, vsExpm1)
+FUNC_FF(FUNC_ABS_FF,     "absolute_ff", fabsf,  fabsf2,  vsAbs)
+FUNC_FF(FUNC_CONJ_FF,    "conjugate_ff",fconjf, fconjf2, vsConj)
+FUNC_FF(FUNC_CEIL_FF,    "ceil_ff",     ceilf,  ceilf2,  vsCeil)
+FUNC_FF(FUNC_FLOOR_FF,   "floor_ff",    floorf, floorf2, vsFloor)
+FUNC_FF(FUNC_TRUNC_FF,   "trunc_ff",    truncf, truncf2, vsTrunc)
+FUNC_FF(FUNC_SIGN_FF,   "sign_ff",    signf, signf2, vsSign)
+//rint rounds to nearest even integer, matching NumPy (round doesn't)
+FUNC_FF(FUNC_ROUND_FF,   "round_ff",    rintf,  rintf2,  vsRint)
+FUNC_FF(FUNC_FF_LAST,    NULL,          NULL,   NULL,    NULL)
+#ifdef ELIDE_FUNC_FF
+#undef ELIDE_FUNC_FF
+#undef FUNC_FF
+#endif
+
+#ifndef FUNC_FFF
+#define ELIDE_FUNC_FFF
+#define FUNC_FFF(...)
+#endif
+FUNC_FFF(FUNC_FMOD_FFF,    "fmod_fff",    fmodf,  fmodf2,  vsfmod)
+FUNC_FFF(FUNC_ARCTAN2_FFF, "arctan2_fff", atan2f, atan2f2, vsAtan2)
+FUNC_FFF(FUNC_HYPOT_FFF, "hypot_fff", hypotf, hypotf2, vsHypot)
+FUNC_FFF(FUNC_NEXTAFTER_FFF, "nextafter_fff", nextafterf, nextafterf2, vsNextAfter)
+FUNC_FFF(FUNC_COPYSIGN_FFF, "copysign_fff", copysignf, copysignf2, vsCopySign)
+FUNC_FFF(FUNC_MAXIMUM_FFF,   "maximum_fff",    fmaxf_, fmaxf2, vsFmax_)
+FUNC_FFF(FUNC_MINIMUM_FFF,   "minimum_fff",    fminf_, fminf2, vsFmin_)
+FUNC_FFF(FUNC_FFF_LAST,    NULL,          NULL,   NULL,    NULL)
+#ifdef ELIDE_FUNC_FFF
+#undef ELIDE_FUNC_FFF
+#undef FUNC_FFF
+#endif
+
+#ifndef FUNC_DD
+#define ELIDE_FUNC_DD
+#define FUNC_DD(...)
+#endif
+FUNC_DD(FUNC_SQRT_DD,    "sqrt_dd",     sqrt,  vdSqrt)
+FUNC_DD(FUNC_SIN_DD,     "sin_dd",      sin,   vdSin)
+FUNC_DD(FUNC_COS_DD,     "cos_dd",      cos,   vdCos)
+FUNC_DD(FUNC_TAN_DD,     "tan_dd",      tan,   vdTan)
+FUNC_DD(FUNC_ARCSIN_DD,  "arcsin_dd",   asin,  vdAsin)
+FUNC_DD(FUNC_ARCCOS_DD,  "arccos_dd",   acos,  vdAcos)
+FUNC_DD(FUNC_ARCTAN_DD,  "arctan_dd",   atan,  vdAtan)
+FUNC_DD(FUNC_SINH_DD,    "sinh_dd",     sinh,  vdSinh)
+FUNC_DD(FUNC_COSH_DD,    "cosh_dd",     cosh,  vdCosh)
+FUNC_DD(FUNC_TANH_DD,    "tanh_dd",     tanh,  vdTanh)
+FUNC_DD(FUNC_ARCSINH_DD, "arcsinh_dd",  asinh, vdAsinh)
+FUNC_DD(FUNC_ARCCOSH_DD, "arccosh_dd",  acosh, vdAcosh)
+FUNC_DD(FUNC_ARCTANH_DD, "arctanh_dd",  atanh, vdAtanh)
+FUNC_DD(FUNC_LOG_DD,     "log_dd",      log,   vdLn)
+FUNC_DD(FUNC_LOG1P_DD,   "log1p_dd",    log1p, vdLog1p)
+FUNC_DD(FUNC_LOG10_DD,   "log10_dd",    log10, vdLog10)
+FUNC_DD(FUNC_LOG2_DD,   "log2_dd",      log2,  vdLog2)
+FUNC_DD(FUNC_EXP_DD,     "exp_dd",      exp,   vdExp)
+FUNC_DD(FUNC_EXPM1_DD,   "expm1_dd",    expm1, vdExpm1)
+FUNC_DD(FUNC_ABS_DD,     "absolute_dd", fabs,  vdAbs)
+FUNC_DD(FUNC_CONJ_DD,    "conjugate_dd",fconj, vdConj)
+FUNC_DD(FUNC_CEIL_DD,    "ceil_dd",     ceil,  vdCeil)
+FUNC_DD(FUNC_FLOOR_DD,   "floor_dd",    floor, vdFloor)
+FUNC_DD(FUNC_TRUNC_DD,   "trunc_dd",    trunc, vdTrunc)
+FUNC_DD(FUNC_SIGN_DD,    "sign_dd",     sign,  vdSign)
+//rint rounds to nearest even integer, matching NumPy (round doesn't)
+FUNC_DD(FUNC_ROUND_DD,   "round_dd",    rint,  vdRint)
+FUNC_DD(FUNC_DD_LAST,    NULL,          NULL,  NULL)
+#ifdef ELIDE_FUNC_DD
+#undef ELIDE_FUNC_DD
+#undef FUNC_DD
+#endif
+
+// double -> boolean functions
+#ifndef FUNC_BD
+#define ELIDE_FUNC_BD
+#define FUNC_BD(...)
+#endif
+FUNC_BD(FUNC_ISNAN_BD,   "isnan_bd",    isnand, vdIsnan)
+FUNC_BD(FUNC_ISFINITE_BD, "isfinite_bd", isfinited, vdIsfinite)
+FUNC_BD(FUNC_ISINF_BD, "isinf_bd", isinfd, vdIsinf)
+FUNC_BD(FUNC_SIGNBIT_BD, "signbit_bd",  signbit, vdSignBit)
+FUNC_BD(FUNC_BD_LAST,    NULL,          NULL,  NULL)
+#ifdef ELIDE_FUNC_BD
+#undef ELIDE_FUNC_BD
+#undef FUNC_BD
+#endif
+
+// float -> boolean functions (C99 defines the same function for all types)
+#ifndef FUNC_BF
+#define ELIDE_FUNC_BF
+#define FUNC_BF(...)
+#endif // use wrappers as there is name collision with isnanf in std
+FUNC_BF(FUNC_ISNAN_BF,   "isnan_bf",    isnanf_,  isnanf2, vsIsnan)
+FUNC_BF(FUNC_ISFINITE_BF, "isfinite_bf", isfinitef_, isfinitef2, vsIsfinite)
+FUNC_BF(FUNC_ISINF_BF, "isinf_bf", isinff_, isinff2, vsIsinf)
+FUNC_BF(FUNC_SIGNBIT_BF, "signbit_bf", signbitf, signbitf2, vsSignBit)
+FUNC_BF(FUNC_BF_LAST,    NULL,            NULL,        NULL,    NULL)
+#ifdef ELIDE_FUNC_BF
+#undef ELIDE_FUNC_BF
+#undef FUNC_BF
+#endif
+
+#ifndef FUNC_DDD
+#define ELIDE_FUNC_DDD
+#define FUNC_DDD(...)
+#endif
+FUNC_DDD(FUNC_FMOD_DDD,    "fmod_ddd",    fmod,  vdfmod)
+FUNC_DDD(FUNC_ARCTAN2_DDD, "arctan2_ddd", atan2, vdAtan2)
+FUNC_DDD(FUNC_HYPOT_DDD, "hypot_ddd", hypot, vdHypot)
+FUNC_DDD(FUNC_NEXTAFTER_DDD, "nextafter_ddd", nextafter, vdNextAfter)
+FUNC_DDD(FUNC_COPYSIGN_DDD, "copysign_ddd", copysign, vdCopySign)
+FUNC_DDD(FUNC_MAXIMUM_DDD, "maximum_ddd",  fmaxd,  vdFmax_)
+FUNC_DDD(FUNC_MINIMUM_DDD, "minimum_ddd",  fmind,  vdFmin_)
+FUNC_DDD(FUNC_DDD_LAST,    NULL,          NULL,  NULL)
+#ifdef ELIDE_FUNC_DDD
+#undef ELIDE_FUNC_DDD
+#undef FUNC_DDD
+#endif
+
+#ifndef FUNC_CC
+#define ELIDE_FUNC_CC
+#define FUNC_CC(...)
+#endif
+FUNC_CC(FUNC_SQRT_CC,    "sqrt_cc",     nc_sqrt,   vzSqrt)
+FUNC_CC(FUNC_SIN_CC,     "sin_cc",      nc_sin,    vzSin)
+FUNC_CC(FUNC_COS_CC,     "cos_cc",      nc_cos,    vzCos)
+FUNC_CC(FUNC_TAN_CC,     "tan_cc",      nc_tan,    vzTan)
+FUNC_CC(FUNC_ARCSIN_CC,  "arcsin_cc",   nc_asin,   vzAsin)
+FUNC_CC(FUNC_ARCCOS_CC,  "arccos_cc",   nc_acos,   vzAcos)
+FUNC_CC(FUNC_ARCTAN_CC,  "arctan_cc",   nc_atan,   vzAtan)
+FUNC_CC(FUNC_SINH_CC,    "sinh_cc",     nc_sinh,   vzSinh)
+FUNC_CC(FUNC_COSH_CC,    "cosh_cc",     nc_cosh,   vzCosh)
+FUNC_CC(FUNC_TANH_CC,    "tanh_cc",     nc_tanh,   vzTanh)
+FUNC_CC(FUNC_ARCSINH_CC, "arcsinh_cc",  nc_asinh,  vzAsinh)
+FUNC_CC(FUNC_ARCCOSH_CC, "arccosh_cc",  nc_acosh,  vzAcosh)
+FUNC_CC(FUNC_ARCTANH_CC, "arctanh_cc",  nc_atanh,  vzAtanh)
+FUNC_CC(FUNC_LOG_CC,     "log_cc",      nc_log,    vzLn)
+FUNC_CC(FUNC_LOG1P_CC,   "log1p_cc",    nc_log1p,  vzLog1p)
+FUNC_CC(FUNC_LOG10_CC,   "log10_cc",    nc_log10,  vzLog10)
+FUNC_CC(FUNC_LOG2_CC,    "log2_cc",     nc_log2,   vzLog2)
+FUNC_CC(FUNC_EXP_CC,     "exp_cc",      nc_exp,    vzExp)
+FUNC_CC(FUNC_EXPM1_CC,   "expm1_cc",    nc_expm1,  vzExpm1)
+FUNC_CC(FUNC_ABS_CC,     "absolute_cc", nc_abs,    vzAbs_)
+FUNC_CC(FUNC_CONJ_CC,    "conjugate_cc",nc_conj,   vzConj)
+FUNC_CC(FUNC_SIGN_CC,   "sign_cc",    nc_sign,   vzSign)
+// rint rounds to nearest even integer, matches NumPy behaviour (round doesn't)
+FUNC_CC(FUNC_ROUND_CC,   "round_cc",    nc_rint,   vzRint)
+FUNC_CC(FUNC_CC_LAST,    NULL,          NULL,      NULL)
+#ifdef ELIDE_FUNC_CC
+#undef ELIDE_FUNC_CC
+#undef FUNC_CC
+#endif
+
+#ifndef FUNC_CCC
+#define ELIDE_FUNC_CCC
+#define FUNC_CCC(...)
+#endif
+FUNC_CCC(FUNC_POW_CCC,   "pow_ccc", nc_pow)
+FUNC_CCC(FUNC_CCC_LAST,  NULL,      NULL)
+#ifdef ELIDE_FUNC_CCC
+#undef ELIDE_FUNC_CCC
+#undef FUNC_CCC
+#endif
+
+// complex -> boolean functions
+#ifndef FUNC_BC
+#define ELIDE_FUNC_BC
+#define FUNC_BC(...)
+#endif // use wrappers as there is name collision with isnanf in std
+FUNC_BC(FUNC_ISNAN_BC,   "isnan_bc",    nc_isnan, vzIsnan)
+FUNC_BC(FUNC_ISFINITE_BC, "isfinite_bc", nc_isfinite, vzIsfinite)
+FUNC_BC(FUNC_ISINF_BC, "isinf_bc", nc_isinf, vzIsinf)
+FUNC_BC(FUNC_BC_LAST,    NULL,            NULL,        NULL)
+#ifdef ELIDE_FUNC_BC
+#undef ELIDE_FUNC_BC
+#undef FUNC_BC
+#endif
+
+// int -> int functions
+#ifndef FUNC_II
+#define ELIDE_FUNC_II
+#define FUNC_II(...)
+#endif
+FUNC_II(FUNC_SIGN_II,   "sign_ii",  signi, viSign)
+FUNC_II(FUNC_ROUND_II,  "round_ii", rinti, viRint)
+FUNC_II(FUNC_ABS_II,  "absolute_ii", fabsi, viFabs)
+FUNC_II(FUNC_II_LAST,  NULL,      NULL, NULL)
+#ifdef ELIDE_FUNC_II
+#undef ELIDE_FUNC_II
+#undef FUNC_II
+#endif
+
+#ifndef FUNC_LL
+#define ELIDE_FUNC_LL
+#define FUNC_LL(...)
+#endif
+FUNC_LL(FUNC_SIGN_LL,  "sign_ll",   signl,   vlSign)
+FUNC_LL(FUNC_ROUND_LL, "round_ll",  rintl,   vlRint)
+FUNC_LL(FUNC_ABS_LL,  "absolute_ll", fabsl, vlFabs)
+FUNC_LL(FUNC_LL_LAST,  NULL,      NULL, NULL)
+#ifdef ELIDE_FUNC_LL
+#undef ELIDE_FUNC_LL
+#undef FUNC_LL
+#endif

numexpr/interp_body.cpp

@@ -0,0 +1,602 @@
+/*********************************************************************
+  Numexpr - Fast numerical array expression evaluator for NumPy.
+
+      License: MIT
+      Author:  See AUTHORS.txt
+
+  See LICENSE.txt for details about copyright and rights to use.
+**********************************************************************/
+
+// WARNING: This file is included multiple times in `interpreter.cpp`. It is
+// essentially a very macro-heavy jump table. Interpretation is best done by
+// the developer by expanding all macros (e.g. adding `'-E'` to the `extra_cflags`
+// argument in `setup.py` and looking at the resulting `interpreter.cpp`.
+//
+// Changes made to this file will not be recognized by the compile, so the developer
+// must make a trivial change is made to `interpreter.cpp` or delete the `build/`
+// directory in-between each build.
+{
+#define VEC_LOOP(expr) for(j = 0; j < BLOCK_SIZE; j++) {       \
+        expr;                                   \
+    }
+
+#define VEC_ARG0(expr)                          \
+    BOUNDS_CHECK(store_in);                     \
+    {                                           \
+        char *dest = mem[store_in];             \
+        VEC_LOOP(expr);                         \
+    } break
+
+#define VEC_ARG1(expr)                          \
+    BOUNDS_CHECK(store_in);                     \
+    BOUNDS_CHECK(arg1);                         \
+    {                                           \
+        char *dest = mem[store_in];             \
+        char *x1 = mem[arg1];                   \
+        npy_intp ss1 = params.memsizes[arg1];       \
+        npy_intp sb1 = memsteps[arg1];              \
+        /* nowarns is defined and used so as to \
+        avoid compiler warnings about unused    \
+        variables */                            \
+        npy_intp nowarns = ss1+sb1+*x1;             \
+        nowarns += 1;                           \
+        VEC_LOOP(expr);                         \
+    } break
+
+#define VEC_ARG2(expr)                          \
+    BOUNDS_CHECK(store_in);                     \
+    BOUNDS_CHECK(arg1);                         \
+    BOUNDS_CHECK(arg2);                         \
+    {                                           \
+        char *dest = mem[store_in];             \
+        char *x1 = mem[arg1];                   \
+        npy_intp ss1 = params.memsizes[arg1];       \
+        npy_intp sb1 = memsteps[arg1];              \
+        /* nowarns is defined and used so as to \
+        avoid compiler warnings about unused    \
+        variables */                            \
+        npy_intp nowarns = ss1+sb1+*x1;             \
+        char *x2 = mem[arg2];                   \
+        npy_intp ss2 = params.memsizes[arg2];       \
+        npy_intp sb2 = memsteps[arg2];              \
+        nowarns += ss2+sb2+*x2;                 \
+        VEC_LOOP(expr);                         \
+    } break
+
+#define VEC_ARG3(expr)                          \
+    BOUNDS_CHECK(store_in);                     \
+    BOUNDS_CHECK(arg1);                         \
+    BOUNDS_CHECK(arg2);                         \
+    BOUNDS_CHECK(arg3);                         \
+    {                                           \
+        char *dest = mem[store_in];             \
+        char *x1 = mem[arg1];                   \
+        npy_intp ss1 = params.memsizes[arg1];       \
+        npy_intp sb1 = memsteps[arg1];              \
+        /* nowarns is defined and used so as to \
+        avoid compiler warnings about unused    \
+        variables */                            \
+        npy_intp nowarns = ss1+sb1+*x1;             \
+        char *x2 = mem[arg2];                   \
+        npy_intp ss2 = params.memsizes[arg2];       \
+        npy_intp sb2 = memsteps[arg2];              \
+        char *x3 = mem[arg3];                   \
+        npy_intp ss3 = params.memsizes[arg3];       \
+        npy_intp sb3 = memsteps[arg3];              \
+        nowarns += ss2+sb2+*x2;                 \
+        nowarns += ss3+sb3+*x3;                 \
+        VEC_LOOP(expr);                         \
+    } break
+
+#define VEC_ARG1_VML(expr)                      \
+    BOUNDS_CHECK(store_in);                     \
+    BOUNDS_CHECK(arg1);                         \
+    {                                           \
+        char *dest = mem[store_in];             \
+        char *x1 = mem[arg1];                   \
+        expr;                                   \
+    } break
+
+#define VEC_ARG2_VML(expr)                      \
+    BOUNDS_CHECK(store_in);                     \
+    BOUNDS_CHECK(arg1);                         \
+    BOUNDS_CHECK(arg2);                         \
+    {                                           \
+        char *dest = mem[store_in];             \
+        char *x1 = mem[arg1];                   \
+        char *x2 = mem[arg2];                   \
+        expr;                                   \
+    } break
+
+#define VEC_ARG3_VML(expr)                      \
+    BOUNDS_CHECK(store_in);                     \
+    BOUNDS_CHECK(arg1);                         \
+    BOUNDS_CHECK(arg2);                         \
+    BOUNDS_CHECK(arg3);                         \
+    {                                           \
+        char *dest = mem[store_in];             \
+        char *x1 = mem[arg1];                   \
+        char *x2 = mem[arg2];                   \
+        char *x3 = mem[arg3];                   \
+        expr;                                   \
+    } break
+
+    int pc;
+    unsigned int j;
+
+    // set up pointers to next block of inputs and outputs
+#ifdef SINGLE_ITEM_CONST_LOOP
+    mem[0] = params.output;
+#else // SINGLE_ITEM_CONST_LOOP
+    // use the iterator's inner loop data
+    memcpy(mem, iter_dataptr, (1+params.n_inputs)*sizeof(char*));
+#  ifndef NO_OUTPUT_BUFFERING
+    // if output buffering is necessary, first write to the buffer
+    if(params.out_buffer != NULL) {
+        mem[0] = params.out_buffer;
+    }
+#  endif // NO_OUTPUT_BUFFERING
+    memcpy(memsteps, iter_strides, (1+params.n_inputs)*sizeof(npy_intp));
+#endif // SINGLE_ITEM_CONST_LOOP
+
+    // WARNING: From now on, only do references to mem[arg[123]]
+    // & memsteps[arg[123]] inside the VEC_ARG[123] macros,
+    // or you will risk accessing invalid addresses.
+
+    for (pc = 0; pc < params.prog_len; pc += 4) {
+        unsigned char op = params.program[pc];
+        unsigned int store_in = params.program[pc+1];
+        unsigned int arg1 = params.program[pc+2];
+        unsigned int arg2 = params.program[pc+3];
+        #define      arg3   params.program[pc+5]
+        // Iterator reduce macros
+#ifdef REDUCTION_INNER_LOOP // Reduce is the inner loop
+        #define i_reduce    *(int *)dest
+        #define l_reduce    *(long long *)dest
+        #define f_reduce    *(float *)dest
+        #define d_reduce    *(double *)dest
+        #define cr_reduce   *(double *)dest
+        #define ci_reduce   *((double *)dest+1)
+#else /* Reduce is the outer loop */
+        #define i_reduce    i_dest
+        #define l_reduce    l_dest
+        #define f_reduce    f_dest
+        #define d_reduce    d_dest
+        #define cr_reduce   cr_dest
+        #define ci_reduce   ci_dest
+#endif
+        #define b_dest ((char *)dest)[j]
+        #define i_dest ((int *)dest)[j]
+        #define l_dest ((long long *)dest)[j]
+        #define f_dest ((float *)dest)[j]
+        #define d_dest ((double *)dest)[j]
+        #define cr_dest ((double *)dest)[2*j]
+        #define ci_dest ((double *)dest)[2*j+1]
+        #define s_dest ((char *)dest + j*memsteps[store_in])
+        #define b1    ((char   *)(x1+j*sb1))[0]
+        #define i1    ((int    *)(x1+j*sb1))[0]
+        #define l1    ((long long *)(x1+j*sb1))[0]
+        #define f1    ((float  *)(x1+j*sb1))[0]
+        #define d1    ((double *)(x1+j*sb1))[0]
+        #define c1r   ((double *)(x1+j*sb1))[0]
+        #define c1i   ((double *)(x1+j*sb1))[1]
+        #define s1    ((char   *)x1+j*sb1)
+        #define b2    ((char   *)(x2+j*sb2))[0]
+        #define i2    ((int    *)(x2+j*sb2))[0]
+        #define l2    ((long long *)(x2+j*sb2))[0]
+        #define f2    ((float  *)(x2+j*sb2))[0]
+        #define d2    ((double *)(x2+j*sb2))[0]
+        #define c2r   ((double *)(x2+j*sb2))[0]
+        #define c2i   ((double *)(x2+j*sb2))[1]
+        #define s2    ((char   *)x2+j*sb2)
+        #define b3    ((char   *)(x3+j*sb3))[0]
+        #define i3    ((int    *)(x3+j*sb3))[0]
+        #define l3    ((long long *)(x3+j*sb3))[0]
+        #define f3    ((float  *)(x3+j*sb3))[0]
+        #define d3    ((double *)(x3+j*sb3))[0]
+        #define c3r   ((double *)(x3+j*sb3))[0]
+        #define c3i   ((double *)(x3+j*sb3))[1]
+        #define s3    ((char   *)x3+j*sb3)
+        /* Some temporaries */
+        double da, db;
+        std::complex<double> ca, cb;
+
+        switch (op) {
+
+        case OP_NOOP: break;
+
+        case OP_COPY_BB: VEC_ARG1(b_dest = b1);
+        case OP_COPY_SS: VEC_ARG1(memcpy(s_dest, s1, ss1));
+        /* The next versions of copy opcodes can cope with unaligned
+           data even on platforms that crash while accessing it
+           (like the Sparc architecture under Solaris). */
+        case OP_COPY_II: VEC_ARG1(memcpy(&i_dest, s1, sizeof(int)));
+        case OP_COPY_LL: VEC_ARG1(memcpy(&l_dest, s1, sizeof(long long)));
+        case OP_COPY_FF: VEC_ARG1(memcpy(&f_dest, s1, sizeof(float)));
+        case OP_COPY_DD: VEC_ARG1(memcpy(&d_dest, s1, sizeof(double)));
+        case OP_COPY_CC: VEC_ARG1(memcpy(&cr_dest, s1, sizeof(double)*2));
+
+        /* Bool */
+        case OP_INVERT_BB: VEC_ARG1(b_dest = !b1);
+        case OP_AND_BBB: VEC_ARG2(b_dest = (b1 && b2));
+        case OP_OR_BBB: VEC_ARG2(b_dest = (b1 || b2));
+        case OP_XOR_BBB: VEC_ARG2(b_dest = (b1 || b2) && !(b1 && b2) );
+
+        case OP_EQ_BBB: VEC_ARG2(b_dest = (b1 == b2));
+        case OP_NE_BBB: VEC_ARG2(b_dest = (b1 != b2));
+        case OP_WHERE_BBBB: VEC_ARG3(b_dest = b1 ? b2 : b3);
+
+        /* Comparisons */
+        case OP_GT_BII: VEC_ARG2(b_dest = (i1 > i2));
+        case OP_GE_BII: VEC_ARG2(b_dest = (i1 >= i2));
+        case OP_EQ_BII: VEC_ARG2(b_dest = (i1 == i2));
+        case OP_NE_BII: VEC_ARG2(b_dest = (i1 != i2));
+
+        case OP_GT_BLL: VEC_ARG2(b_dest = (l1 > l2));
+        case OP_GE_BLL: VEC_ARG2(b_dest = (l1 >= l2));
+        case OP_EQ_BLL: VEC_ARG2(b_dest = (l1 == l2));
+        case OP_NE_BLL: VEC_ARG2(b_dest = (l1 != l2));
+
+        case OP_GT_BFF: VEC_ARG2(b_dest = (f1 > f2));
+        case OP_GE_BFF: VEC_ARG2(b_dest = (f1 >= f2));
+        case OP_EQ_BFF: VEC_ARG2(b_dest = (f1 == f2));
+        case OP_NE_BFF: VEC_ARG2(b_dest = (f1 != f2));
+
+        case OP_GT_BDD: VEC_ARG2(b_dest = (d1 > d2));
+        case OP_GE_BDD: VEC_ARG2(b_dest = (d1 >= d2));
+        case OP_EQ_BDD: VEC_ARG2(b_dest = (d1 == d2));
+        case OP_NE_BDD: VEC_ARG2(b_dest = (d1 != d2));
+
+        case OP_GT_BSS: VEC_ARG2(b_dest = (stringcmp(s1, s2, ss1, ss2) > 0));
+        case OP_GE_BSS: VEC_ARG2(b_dest = (stringcmp(s1, s2, ss1, ss2) >= 0));
+        case OP_EQ_BSS: VEC_ARG2(b_dest = (stringcmp(s1, s2, ss1, ss2) == 0));
+        case OP_NE_BSS: VEC_ARG2(b_dest = (stringcmp(s1, s2, ss1, ss2) != 0));
+
+        case OP_CONTAINS_BSS: VEC_ARG2(b_dest = stringcontains(s1, s2, ss1, ss2));
+
+        /* Int */
+        case OP_CAST_IB: VEC_ARG1(i_dest = (int)(b1));
+        case OP_ONES_LIKE_II: VEC_ARG0(i_dest = 1);
+        case OP_NEG_II: VEC_ARG1(i_dest = -i1);
+
+        case OP_ADD_III: VEC_ARG2(i_dest = i1 + i2);
+        case OP_SUB_III: VEC_ARG2(i_dest = i1 - i2);
+        case OP_MUL_III: VEC_ARG2(i_dest = i1 * i2);
+        case OP_DIV_III: VEC_ARG2(i_dest = i2 ? (i1 / i2) : 0);
+        case OP_POW_III: VEC_ARG2(i_dest = (i2 < 0) ? (1 / i1) : (int)pow((double)i1, i2));
+        case OP_MOD_III: VEC_ARG2(i_dest = i2 == 0 ? 0 :((i1 % i2) + i2) % i2);
+        case OP_FLOORDIV_III: VEC_ARG2(i_dest = i2 ? (i1 / i2) - ((i1 % i2 != 0) && (i1 < 0 != i2 < 0)) : 0);
+        case OP_LSHIFT_III: VEC_ARG2(i_dest = i1 << i2);
+        case OP_RSHIFT_III: VEC_ARG2(i_dest = i1 >> i2);
+
+        case OP_WHERE_IBII: VEC_ARG3(i_dest = b1 ? i2 : i3);
+        //Bitwise ops
+        case OP_INVERT_II: VEC_ARG1(i_dest = ~i1);
+        case OP_AND_III: VEC_ARG2(i_dest = (i1 & i2));
+        case OP_OR_III: VEC_ARG2(i_dest = (i1 | i2));
+        case OP_XOR_III: VEC_ARG2(i_dest = (i1 ^ i2));
+
+        /* Long */
+        case OP_CAST_LI: VEC_ARG1(l_dest = (long long)(i1));
+        case OP_ONES_LIKE_LL: VEC_ARG0(l_dest = 1);
+        case OP_NEG_LL: VEC_ARG1(l_dest = -l1);
+
+        case OP_ADD_LLL: VEC_ARG2(l_dest = l1 + l2);
+        case OP_SUB_LLL: VEC_ARG2(l_dest = l1 - l2);
+        case OP_MUL_LLL: VEC_ARG2(l_dest = l1 * l2);
+        case OP_DIV_LLL: VEC_ARG2(l_dest = l2 ? (l1 / l2) : 0);
+#if defined _MSC_VER && _MSC_VER < 1800
+        case OP_POW_LLL: VEC_ARG2(l_dest = (l2 < 0) ? (1 / l1) : (long long)pow((long double)l1, (long double)l2));
+#else
+        case OP_POW_LLL: VEC_ARG2(l_dest = (l2 < 0) ? (1 / l1) : (long long)llround(pow((long double)l1, (long double)l2)));
+#endif
+        case OP_MOD_LLL: VEC_ARG2(l_dest = l2 == 0 ? 0 :((l1 % l2) + l2) % l2);
+        case OP_FLOORDIV_LLL: VEC_ARG2(l_dest = l2 ? (l1 / l2) - ((l1 % l2 != 0) && (l1 < 0 != l2 < 0)): 0);
+        case OP_LSHIFT_LLL: VEC_ARG2(l_dest = l1 << l2);
+        case OP_RSHIFT_LLL: VEC_ARG2(l_dest = l1 >> l2);
+
+        case OP_WHERE_LBLL: VEC_ARG3(l_dest = b1 ? l2 : l3);
+        //Bitwise ops
+        case OP_INVERT_LL: VEC_ARG1(l_dest = ~l1);
+        case OP_AND_LLL: VEC_ARG2(l_dest = (l1 & l2));
+        case OP_OR_LLL: VEC_ARG2(l_dest = (l1 | l2));
+        case OP_XOR_LLL: VEC_ARG2(l_dest = (l1 ^ l2));
+
+        /* Float */
+        case OP_CAST_FI: VEC_ARG1(f_dest = (float)(i1));
+        case OP_CAST_FL: VEC_ARG1(f_dest = (float)(l1));
+        case OP_ONES_LIKE_FF: VEC_ARG0(f_dest = 1.0);
+        case OP_NEG_FF: VEC_ARG1(f_dest = -f1);
+
+        case OP_ADD_FFF: VEC_ARG2(f_dest = f1 + f2);
+        case OP_SUB_FFF: VEC_ARG2(f_dest = f1 - f2);
+        case OP_MUL_FFF: VEC_ARG2(f_dest = f1 * f2);
+        case OP_DIV_FFF:
+#ifdef USE_VML
+            VEC_ARG2_VML(vsDiv(BLOCK_SIZE,
+                               (float*)x1, (float*)x2, (float*)dest));
+#else
+            VEC_ARG2(f_dest = f1 / f2);
+#endif
+        case OP_POW_FFF:
+#ifdef USE_VML
+            VEC_ARG2_VML(vsPow(BLOCK_SIZE,
+                               (float*)x1, (float*)x2, (float*)dest));
+#else
+            VEC_ARG2(f_dest = powf(f1, f2));
+#endif
+        case OP_MOD_FFF: VEC_ARG2(f_dest = f1 - floorf(f1/f2) * f2);
+        case OP_FLOORDIV_FFF: VEC_ARG2(f_dest = floorf(f1/f2));
+
+        case OP_SQRT_FF:
+#ifdef USE_VML
+            VEC_ARG1_VML(vsSqrt(BLOCK_SIZE, (float*)x1, (float*)dest));
+#else
+            VEC_ARG1(f_dest = sqrtf(f1));
+#endif
+
+        case OP_WHERE_FBFF: VEC_ARG3(f_dest = b1 ? f2 : f3);
+
+        case OP_FUNC_FFN:
+#ifdef USE_VML
+            VEC_ARG1_VML(functions_ff_vml[arg2](BLOCK_SIZE,
+                                                (float*)x1, (float*)dest));
+#else
+            VEC_ARG1(f_dest = functions_ff[arg2](f1));
+#endif
+        case OP_FUNC_FFFN:
+#ifdef USE_VML
+            VEC_ARG2_VML(functions_fff_vml[arg3](BLOCK_SIZE,
+                                                 (float*)x1, (float*)x2,
+                                                 (float*)dest));
+#else
+            VEC_ARG2(f_dest = functions_fff[arg3](f1, f2));
+#endif
+
+        /* Double */
+        case OP_CAST_DI: VEC_ARG1(d_dest = (double)(i1));
+        case OP_CAST_DL: VEC_ARG1(d_dest = (double)(l1));
+        case OP_CAST_DF: VEC_ARG1(d_dest = (double)(f1));
+        case OP_ONES_LIKE_DD: VEC_ARG0(d_dest = 1.0);
+        case OP_NEG_DD: VEC_ARG1(d_dest = -d1);
+
+        case OP_ADD_DDD: VEC_ARG2(d_dest = d1 + d2);
+        case OP_SUB_DDD: VEC_ARG2(d_dest = d1 - d2);
+        case OP_MUL_DDD: VEC_ARG2(d_dest = d1 * d2);
+        case OP_DIV_DDD:
+#ifdef USE_VML
+            VEC_ARG2_VML(vdDiv(BLOCK_SIZE,
+                               (double*)x1, (double*)x2, (double*)dest));
+#else
+            VEC_ARG2(d_dest = d1 / d2);
+#endif
+        case OP_POW_DDD:
+#ifdef USE_VML
+            VEC_ARG2_VML(vdPow(BLOCK_SIZE,
+                               (double*)x1, (double*)x2, (double*)dest));
+#else
+            VEC_ARG2(d_dest = pow(d1, d2));
+#endif
+        case OP_MOD_DDD: VEC_ARG2(d_dest = d1 - floor(d1/d2) * d2);
+        case OP_FLOORDIV_DDD: VEC_ARG2(d_dest = floor(d1/d2));
+
+        case OP_SQRT_DD:
+#ifdef USE_VML
+            VEC_ARG1_VML(vdSqrt(BLOCK_SIZE, (double*)x1, (double*)dest));
+#else
+            VEC_ARG1(d_dest = sqrt(d1));
+#endif
+
+        case OP_WHERE_DBDD: VEC_ARG3(d_dest = b1 ? d2 : d3);
+
+        case OP_FUNC_DDN:
+#ifdef USE_VML
+            VEC_ARG1_VML(functions_dd_vml[arg2](BLOCK_SIZE,
+                                                (double*)x1, (double*)dest));
+#else
+            VEC_ARG1(d_dest = functions_dd[arg2](d1));
+#endif
+        case OP_FUNC_DDDN:
+#ifdef USE_VML
+            VEC_ARG2_VML(functions_ddd_vml[arg3](BLOCK_SIZE,
+                                                 (double*)x1, (double*)x2,
+                                                 (double*)dest));
+#else
+            VEC_ARG2(d_dest = functions_ddd[arg3](d1, d2));
+#endif
+
+        /* Complex */
+        case OP_CAST_CI: VEC_ARG1(cr_dest = (double)(i1);
+                                  ci_dest = 0);
+        case OP_CAST_CL: VEC_ARG1(cr_dest = (double)(l1);
+                                  ci_dest = 0);
+        case OP_CAST_CF: VEC_ARG1(cr_dest = f1;
+                                  ci_dest = 0);
+        case OP_CAST_CD: VEC_ARG1(cr_dest = d1;
+                                  ci_dest = 0);
+        case OP_ONES_LIKE_CC: VEC_ARG0(cr_dest = 1;
+                                       ci_dest = 0);
+        case OP_NEG_CC: VEC_ARG1(cr_dest = -c1r;
+                                 ci_dest = -c1i);
+
+        case OP_ADD_CCC: VEC_ARG2(cr_dest = c1r + c2r;
+                                  ci_dest = c1i + c2i);
+        case OP_SUB_CCC: VEC_ARG2(cr_dest = c1r - c2r;
+                                  ci_dest = c1i - c2i);
+        case OP_MUL_CCC: VEC_ARG2(da = c1r*c2r - c1i*c2i;
+                                  ci_dest = c1r*c2i + c1i*c2r;
+                                  cr_dest = da);
+        case OP_DIV_CCC:
+#ifdef USE_VMLXXX /* VML complex division is slower */
+            VEC_ARG2_VML(vzDiv(BLOCK_SIZE, (const MKL_Complex16*)x1,
+                               (const MKL_Complex16*)x2, (MKL_Complex16*)dest));
+#else
+            VEC_ARG2(da = c2r*c2r + c2i*c2i;
+                     db = (c1r*c2r + c1i*c2i) / da;
+                     ci_dest = (c1i*c2r - c1r*c2i) / da;
+                     cr_dest = db);
+#endif
+        case OP_EQ_BCC: VEC_ARG2(b_dest = (c1r == c2r && c1i == c2i));
+        case OP_NE_BCC: VEC_ARG2(b_dest = (c1r != c2r || c1i != c2i));
+
+        case OP_WHERE_CBCC: VEC_ARG3(cr_dest = b1 ? c2r : c3r;
+                                     ci_dest = b1 ? c2i : c3i);
+        case OP_FUNC_CCN:
+#ifdef USE_VML
+            VEC_ARG1_VML(functions_cc_vml[arg2](BLOCK_SIZE,
+                                                (const MKL_Complex16*)x1,
+                                                (MKL_Complex16*)dest));
+#else
+            VEC_ARG1(ca.real(c1r);
+                     ca.imag(c1i);
+                     functions_cc[arg2](&ca, &ca);
+                     cr_dest = ca.real();
+                     ci_dest = ca.imag());
+#endif
+        case OP_FUNC_CCCN: VEC_ARG2(ca.real(c1r);
+                                    ca.imag(c1i);
+                                    cb.real(c2r);
+                                    cb.imag(c2i);
+                                    functions_ccc[arg3](&ca, &cb, &ca);
+                                    cr_dest = ca.real();
+                                    ci_dest = ca.imag());
+
+        case OP_REAL_DC: VEC_ARG1(d_dest = c1r);
+        case OP_IMAG_DC: VEC_ARG1(d_dest = c1i);
+        case OP_COMPLEX_CDD: VEC_ARG2(cr_dest = d1;
+                                      ci_dest = d2);
+
+        // Boolean return types
+        case OP_FUNC_BFN:
+#ifdef USE_VML
+            VEC_ARG1_VML(functions_bf_vml[arg2](BLOCK_SIZE,
+                                                (float*)x1, (bool*)dest));
+#else
+            VEC_ARG1(b_dest = functions_bf[arg2](f1));
+#endif
+
+
+        case OP_FUNC_BDN:
+#ifdef USE_VML
+            VEC_ARG1_VML(functions_bd_vml[arg2](BLOCK_SIZE,
+                                                (double*)x1, (bool*)dest));
+#else
+            VEC_ARG1(b_dest = functions_bd[arg2](d1));
+#endif
+
+        case OP_FUNC_BCN:
+#ifdef USE_VML
+            VEC_ARG1_VML(functions_bc_vml[arg2](BLOCK_SIZE,
+                        (const MKL_Complex16*)x1, (bool*)dest));
+#else
+            VEC_ARG1(ca.real(c1r);
+                     ca.imag(c1i);
+                     b_dest = functions_bc[arg2](&ca));
+#endif
+
+        /* Integer return types */
+         case OP_FUNC_IIN:
+#ifdef USE_VML
+            VEC_ARG1_VML(functions_ii_vml[arg2](BLOCK_SIZE,
+                                                (int*)x1, (int*)dest));
+#else
+            VEC_ARG1(i_dest = functions_ii[arg2](i1));
+#endif
+         case OP_FUNC_LLN:
+#ifdef USE_VML
+            VEC_ARG1_VML(functions_ll_vml[arg2](BLOCK_SIZE,
+                                                (long*)x1, (long*)dest));
+#else
+            VEC_ARG1(l_dest = functions_ll[arg2](l1));
+#endif
+
+        /* Reductions */
+        case OP_SUM_IIN: VEC_ARG1(i_reduce += i1);
+        case OP_SUM_LLN: VEC_ARG1(l_reduce += l1);
+        case OP_SUM_FFN: VEC_ARG1(f_reduce += f1);
+        case OP_SUM_DDN: VEC_ARG1(d_reduce += d1);
+        case OP_SUM_CCN: VEC_ARG1(cr_reduce += c1r;
+                                  ci_reduce += c1i);
+
+        case OP_PROD_IIN: VEC_ARG1(i_reduce *= i1);
+        case OP_PROD_LLN: VEC_ARG1(l_reduce *= l1);
+        case OP_PROD_FFN: VEC_ARG1(f_reduce *= f1);
+        case OP_PROD_DDN: VEC_ARG1(d_reduce *= d1);
+        case OP_PROD_CCN: VEC_ARG1(da = cr_reduce*c1r - ci_reduce*c1i;
+                                   ci_reduce = cr_reduce*c1i + ci_reduce*c1r;
+                                   cr_reduce = da);
+
+        case OP_MIN_IIN: VEC_ARG1(i_reduce = fmin(i_reduce, i1));
+        case OP_MIN_LLN: VEC_ARG1(l_reduce = fmin(l_reduce, l1));
+        case OP_MIN_FFN: VEC_ARG1(f_reduce = fmin(f_reduce, f1));
+        case OP_MIN_DDN: VEC_ARG1(d_reduce = fmin(d_reduce, d1));
+
+        case OP_MAX_IIN: VEC_ARG1(i_reduce = fmax(i_reduce, i1));
+        case OP_MAX_LLN: VEC_ARG1(l_reduce = fmax(l_reduce, l1));
+        case OP_MAX_FFN: VEC_ARG1(f_reduce = fmax(f_reduce, f1));
+        case OP_MAX_DDN: VEC_ARG1(d_reduce = fmax(d_reduce, d1));
+
+        default:
+            *pc_error = pc;
+            return -3;
+            break;
+        }
+    }
+
+#ifndef NO_OUTPUT_BUFFERING
+    // If output buffering was necessary, copy the buffer to the output
+    if(params.out_buffer != NULL) {
+        memcpy(iter_dataptr[0], params.out_buffer, params.memsizes[0] * BLOCK_SIZE);
+    }
+#endif // NO_OUTPUT_BUFFERING
+
+#undef VEC_LOOP
+#undef VEC_ARG1
+#undef VEC_ARG2
+#undef VEC_ARG3
+
+#undef i_reduce
+#undef l_reduce
+#undef f_reduce
+#undef d_reduce
+#undef cr_reduce
+#undef ci_reduce
+#undef b_dest
+#undef i_dest
+#undef l_dest
+#undef f_dest
+#undef d_dest
+#undef cr_dest
+#undef ci_dest
+#undef s_dest
+#undef b1
+#undef i1
+#undef l1
+#undef f1
+#undef d1
+#undef c1r
+#undef c1i
+#undef s1
+#undef b2
+#undef i2
+#undef l2
+#undef f2
+#undef d2
+#undef c2r
+#undef c2i
+#undef s2
+#undef b3
+#undef i3
+#undef l3
+#undef f3
+#undef d3
+#undef c3r
+#undef c3i
+#undef s3
+}
+
+/*
+Local Variables:
+   c-basic-offset: 4
+End:
+*/

numexpr/interpreter.hpp

@@ -0,0 +1,137 @@
+#ifndef NUMEXPR_INTERPRETER_HPP
+#define NUMEXPR_INTERPRETER_HPP
+
+#include "numexpr_config.hpp"
+
+// Forward declaration
+struct NumExprObject;
+
+enum OpCodes {
+#define OPCODE(n, e, ...) e = n,
+#include "opcodes.hpp"
+#undef OPCODE
+};
+
+enum FuncFFCodes {
+#define FUNC_FF(fop, ...) fop,
+#include "functions.hpp"
+#undef FUNC_FF
+};
+
+enum FuncBFCodes {
+#define FUNC_BF(fop, ...) fop,
+#include "functions.hpp"
+#undef FUNC_BF
+};
+
+enum FuncFFFCodes {
+#define FUNC_FFF(fop, ...) fop,
+#include "functions.hpp"
+#undef FUNC_FFF
+};
+
+enum FuncDDCodes {
+#define FUNC_DD(fop, ...) fop,
+#include "functions.hpp"
+#undef FUNC_DD
+};
+
+enum FuncBDCodes {
+#define FUNC_BD(fop, ...) fop,
+#include "functions.hpp"
+#undef FUNC_BD
+};
+
+enum FuncBCCodes {
+#define FUNC_BC(fop, ...) fop,
+#include "functions.hpp"
+#undef FUNC_BC
+};
+
+enum FuncIICodes {
+#define FUNC_II(fop, ...) fop,
+#include "functions.hpp"
+#undef FUNC_II
+};
+
+enum FuncLLCodes {
+#define FUNC_LL(fop, ...) fop,
+#include "functions.hpp"
+#undef FUNC_LL
+};
+
+enum FuncDDDCodes {
+#define FUNC_DDD(fop, ...) fop,
+#include "functions.hpp"
+#undef FUNC_DDD
+};
+
+enum FuncCCCodes {
+#define FUNC_CC(fop, ...) fop,
+#include "functions.hpp"
+#undef FUNC_CC
+};
+
+enum FuncCCCCodes {
+#define FUNC_CCC(fop, ...) fop,
+#include "functions.hpp"
+#undef FUNC_CCC
+};
+
+struct vm_params {
+    int prog_len;
+    unsigned char *program;
+    int n_inputs;
+    int n_constants;
+    int n_temps;
+    unsigned int r_end;
+    char *output;
+    char **inputs;
+    char **mem;
+    npy_intp *memsteps;
+    npy_intp *memsizes;
+    struct index_data *index_data;
+    // Memory for output buffering. If output buffering is unneeded,
+    // it contains NULL.
+    char *out_buffer;
+};
+
+// Structure for parameters in worker threads
+struct thread_data {
+    npy_intp start;
+    npy_intp vlen;
+    npy_intp block_size;
+    vm_params params;
+    int ret_code;
+    int *pc_error;
+    char **errmsg;
+    // NOTE: memsteps, iter, and reduce_iter are arrays, they MUST be allocated
+    // to length `global_max_threads` before module load.
+    // One memsteps array per thread
+    // npy_intp *memsteps[MAX_THREADS];
+    npy_intp **memsteps;
+    // One iterator per thread */
+    // NpyIter *iter[MAX_THREADS];
+    NpyIter **iter;
+    // When doing nested iteration for a reduction
+    // NpyIter *reduce_iter[MAX_THREADS]
+    NpyIter **reduce_iter;
+    // Flag indicating reduction is the outer loop instead of the inner
+    bool reduction_outer_loop;
+    // Flag indicating whether output buffering is needed
+    bool need_output_buffering;
+};
+
+// Global state which holds thread parameters
+extern thread_data th_params;
+
+PyObject *NumExpr_run(NumExprObject *self, PyObject *args, PyObject *kwds);
+
+char get_return_sig(PyObject* program);
+int check_program(NumExprObject *self);
+int get_temps_space(const vm_params& params, char **mem, size_t block_size);
+void free_temps_space(const vm_params& params, char **mem);
+int vm_engine_iter_task(NpyIter *iter, npy_intp *memsteps,
+                    const vm_params& params, int *pc_error, char **errmsg);
+
+#endif // NUMEXPR_INTERPRETER_HPP

numexpr/missing_posix_functions.hpp

@@ -0,0 +1,102 @@
+#ifndef NUMEXPR_MISSING_POSIX_FUNCTIONS_HPP
+#define NUMEXPR_MISSING_POSIX_FUNCTIONS_HPP
+
+/*********************************************************************
+  Numexpr - Fast numerical array expression evaluator for NumPy.
+
+      License: MIT
+      Author:  See AUTHORS.txt
+
+  See LICENSE.txt for details about copyright and rights to use.
+**********************************************************************/
+
+/* These functions are not included in some non-POSIX compilers,
+  like MSVC 7.1  */
+
+
+/* Double precision versions */
+
+inline double log1p(double x)
+{
+    double u = 1.0 + x;
+    if (u == 1.0) {
+        return x;
+    } else {
+        return log(u) * x / (u-1.0);
+    }
+}
+
+inline double expm1(double x)
+{
+    double u = exp(x);
+    if (u == 1.0) {
+        return x;
+    } else if (u-1.0 == -1.0) {
+        return -1;
+    } else {
+        return (u-1.0) * x/log(u);
+    }
+}
+
+inline double asinh(double xx)
+{
+    double x, d;
+    int sign;
+    if (xx < 0.0) {
+        sign = -1;
+        x = -xx;
+    }
+    else {
+        sign = 1;
+        x = xx;
+    }
+    if (x > 1e8) {
+        d = x;
+    } else {
+        d = sqrt(x*x + 1.0);
+    }
+    return sign*log1p(x*(1.0 + x/(d+1.0)));
+}
+
+inline double acosh(double x)
+{
+    return 2*log(sqrt((x+1.0)/2)+sqrt((x-1.0)/2));
+}
+
+inline double atanh(double x)
+{
+    /* This definition is different from that in NumPy 1.3 and follows
+    the convention of MatLab.  This will allow for double checking both
+    approaches. */
+    return 0.5*log((1.0+x)/(1.0-x));
+}
+
+
+/* Single precision versions */
+
+inline float log1pf(float x)
+{
+    return (float) log1p((double)x);
+}
+
+inline float expm1f(float x)
+{
+    return (float) expm1((double)x);
+}
+
+inline float asinhf(float x)
+{
+    return (float) asinh((double)x);
+}
+
+inline float acoshf(float x)
+{
+    return (float) acosh((double)x);
+}
+
+inline float atanhf(float x)
+{
+    return (float) atanh((double)x);
+}
+
+#endif // NUMEXPR_MISSING_POSIX_FUNCTIONS_HPP

numexpr/module.cpp

@@ -0,0 +1,552 @@
+// Numexpr - Fast numerical array expression evaluator for NumPy.
+//
+//      License: MIT
+//      Author:  See AUTHORS.txt
+//
+//  See LICENSE.txt for details about copyright and rights to use.
+//
+// module.cpp contains the CPython-specific module exposure.
+
+#define DO_NUMPY_IMPORT_ARRAY
+
+#include "module.hpp"
+#include <structmember.h>
+#include <vector>
+
+#include <signal.h>
+
+#include "interpreter.hpp"
+#include "numexpr_object.hpp"
+
+using namespace std;
+
+// Global state. The file interpreter.hpp also has some global state
+// in its 'th_params' variable
+global_state gs;
+long global_max_threads=DEFAULT_MAX_THREADS;
+
+/* Do the worker job for a certain thread */
+void *th_worker(void *tidptr)
+{
+    int tid = *(int *)tidptr;
+    /* Parameters for threads */
+    npy_intp start;
+    npy_intp vlen;
+    npy_intp block_size;
+    NpyIter *iter;
+    vm_params params;
+    int *pc_error;
+    int ret;
+    int n_inputs;
+    int n_constants;
+    int n_temps;
+    size_t memsize;
+    char **mem;
+    npy_intp *memsteps;
+    npy_intp istart, iend;
+    char **errmsg;
+    // For output buffering if needed
+    vector<char> out_buffer;
+
+    while (1) {
+
+        /* Sentinels have to be initialised yet */
+        if (tid == 0) {
+            gs.init_sentinels_done = 0;
+        }
+
+        /* Meeting point for all threads (wait for initialization) */
+        pthread_mutex_lock(&gs.count_threads_mutex);
+        if (gs.count_threads < gs.nthreads) {
+            gs.count_threads++;
+            /* Beware of spurious wakeups. See issue pydata/numexpr#306. */
+            do {
+                pthread_cond_wait(&gs.count_threads_cv,
+                                  &gs.count_threads_mutex);
+            } while (!gs.barrier_passed);
+        }
+        else {
+            gs.barrier_passed = 1;
+            pthread_cond_broadcast(&gs.count_threads_cv);
+        }
+        pthread_mutex_unlock(&gs.count_threads_mutex);
+
+        /* Check if thread has been asked to return */
+        if (gs.end_threads) {
+            return(0);
+        }
+
+        /* Get parameters for this thread before entering the main loop */
+        start = th_params.start;
+        vlen = th_params.vlen;
+        block_size = th_params.block_size;
+        params = th_params.params;
+        pc_error = th_params.pc_error;
+
+        // If output buffering is needed, allocate it
+        if (th_params.need_output_buffering) {
+            out_buffer.resize(params.memsizes[0] * BLOCK_SIZE1);
+            params.out_buffer = &out_buffer[0];
+        } else {
+            params.out_buffer = NULL;
+        }
+
+        /* Populate private data for each thread */
+        n_inputs = params.n_inputs;
+        n_constants = params.n_constants;
+        n_temps = params.n_temps;
+        memsize = (1+n_inputs+n_constants+n_temps) * sizeof(char *);
+        /* XXX malloc seems thread safe for POSIX, but for Win? */
+        mem = (char **)malloc(memsize);
+        memcpy(mem, params.mem, memsize);
+
+        errmsg = th_params.errmsg;
+
+        params.mem = mem;
+
+        /* Loop over blocks */
+        pthread_mutex_lock(&gs.count_mutex);
+        if (!gs.init_sentinels_done) {
+            /* Set sentinels and other global variables */
+            gs.gindex = start;
+            istart = gs.gindex;
+            iend = istart + block_size;
+            if (iend > vlen) {
+                iend = vlen;
+            }
+            gs.init_sentinels_done = 1;  /* sentinels have been initialised */
+            gs.giveup = 0;            /* no giveup initially */
+        } else {
+            gs.gindex += block_size;
+            istart = gs.gindex;
+            iend = istart + block_size;
+            if (iend > vlen) {
+                iend = vlen;
+            }
+        }
+        /* Grab one of the iterators */
+        iter = th_params.iter[tid];
+        if (iter == NULL) {
+            th_params.ret_code = -1;
+            gs.giveup = 1;
+        }
+        memsteps = th_params.memsteps[tid];
+        /* Get temporary space for each thread */
+        ret = get_temps_space(params, mem, BLOCK_SIZE1);
+        if (ret < 0) {
+            /* Propagate error to main thread */
+            th_params.ret_code = ret;
+            gs.giveup = 1;
+        }
+        pthread_mutex_unlock(&gs.count_mutex);
+
+        while (istart < vlen && !gs.giveup) {
+            /* Reset the iterator to the range for this task */
+            ret = NpyIter_ResetToIterIndexRange(iter, istart, iend,
+                                                errmsg);
+            /* Execute the task */
+            if (ret >= 0) {
+                ret = vm_engine_iter_task(iter, memsteps, params, pc_error, errmsg);
+            }
+
+            if (ret < 0) {
+                pthread_mutex_lock(&gs.count_mutex);
+                gs.giveup = 1;
+                /* Propagate error to main thread */
+                th_params.ret_code = ret;
+                pthread_mutex_unlock(&gs.count_mutex);
+                break;
+            }
+
+            pthread_mutex_lock(&gs.count_mutex);
+            gs.gindex += block_size;
+            istart = gs.gindex;
+            iend = istart + block_size;
+            if (iend > vlen) {
+                iend = vlen;
+            }
+            pthread_mutex_unlock(&gs.count_mutex);
+        }
+
+        /* Meeting point for all threads (wait for finalization) */
+        pthread_mutex_lock(&gs.count_threads_mutex);
+        if (gs.count_threads > 0) {
+            gs.count_threads--;
+            do {
+                pthread_cond_wait(&gs.count_threads_cv,
+                                  &gs.count_threads_mutex);
+            } while (gs.barrier_passed);
+        }
+        else {
+            gs.barrier_passed = 0;
+            pthread_cond_broadcast(&gs.count_threads_cv);
+        }
+        pthread_mutex_unlock(&gs.count_threads_mutex);
+
+        /* Release resources */
+        free_temps_space(params, mem);
+        free(mem);
+
+    }  /* closes while(1) */
+
+    /* This should never be reached, but anyway */
+    return(0);
+}
+
+/* Initialize threads */
+int init_threads(void)
+{
+    int tid, rc;
+
+    if ( !(gs.nthreads > 1 && (!gs.init_threads_done || gs.pid != getpid())) ) {
+        /* Thread pool must always be initialized once and once only. */
+        return(0);
+    }
+
+    /* Initialize mutex and condition variable objects */
+    pthread_mutex_init(&gs.count_mutex, NULL);
+    pthread_mutex_init(&gs.parallel_mutex, NULL);
+
+    /* Barrier initialization */
+    pthread_mutex_init(&gs.count_threads_mutex, NULL);
+    pthread_cond_init(&gs.count_threads_cv, NULL);
+    gs.count_threads = 0;      /* Reset threads counter */
+    gs.barrier_passed = 0;
+
+    /*
+     * Our worker threads should not deal with signals from the rest of the
+     * application - mask everything temporarily in this thread, so our workers
+     * can inherit that mask
+     */
+    sigset_t sigset_block_all, sigset_restore;
+    rc = sigfillset(&sigset_block_all);
+    if (rc != 0) {
+        fprintf(stderr, "ERROR; failed to block signals: sigfillset: %s",
+                strerror(rc));
+        exit(-1);
+    }
+    rc = pthread_sigmask( SIG_BLOCK, &sigset_block_all, &sigset_restore);
+    if (rc != 0) {
+        fprintf(stderr, "ERROR; failed to block signals: pthread_sigmask: %s",
+                strerror(rc));
+        exit(-1);
+    }
+
+    /* Now create the threads */
+    for (tid = 0; tid < gs.nthreads; tid++) {
+        gs.tids[tid] = tid;
+        rc = pthread_create(&gs.threads[tid], NULL, th_worker,
+                            (void *)&gs.tids[tid]);
+        if (rc) {
+            fprintf(stderr,
+                    "ERROR; return code from pthread_create() is %d\n", rc);
+            fprintf(stderr, "\tError detail: %s\n", strerror(rc));
+            exit(-1);
+        }
+    }
+
+    /*
+     * Restore the signal mask so the main thread can process signals as
+     * expected
+     */
+    rc = pthread_sigmask( SIG_SETMASK, &sigset_restore, NULL);
+    if (rc != 0) {
+        fprintf(stderr,
+                "ERROR: failed to restore signal mask: pthread_sigmask: %s",
+                strerror(rc));
+        exit(-1);
+    }
+
+    gs.init_threads_done = 1;                 /* Initialization done! */
+    gs.pid = (int)getpid();                   /* save the PID for this process */
+
+    return(0);
+}
+
+/* Set the number of threads in numexpr's VM */
+int numexpr_set_nthreads(int nthreads_new)
+{
+    int nthreads_old = gs.nthreads;
+    int t, rc;
+    void *status;
+
+    // if (nthreads_new > MAX_THREADS) {
+    //     fprintf(stderr,
+    //             "Error.  nthreads cannot be larger than MAX_THREADS (%d)",
+    //             MAX_THREADS);
+    //     return -1;
+    // }
+    if (nthreads_new > global_max_threads) {
+        fprintf(stderr,
+                "Error.  nthreads cannot be larger than environment variable \"NUMEXPR_MAX_THREADS\" (%ld)",
+                global_max_threads);
+        return -1;
+    }
+    else if (nthreads_new <= 0) {
+        fprintf(stderr, "Error.  nthreads must be a positive integer");
+        return -1;
+    }
+
+    /* Only join threads if they are not initialized or if our PID is
+       different from that in pid var (probably means that we are a
+       subprocess, and thus threads are non-existent). */
+    if (gs.nthreads > 1 && gs.init_threads_done && gs.pid == getpid()) {
+        /* Tell all existing threads to finish */
+        gs.end_threads = 1;
+        pthread_mutex_lock(&gs.count_threads_mutex);
+        if (gs.count_threads < gs.nthreads) {
+            gs.count_threads++;
+            do {
+                pthread_cond_wait(&gs.count_threads_cv,
+                                  &gs.count_threads_mutex);
+            } while (!gs.barrier_passed);
+        }
+        else {
+            gs.barrier_passed = 1;
+            pthread_cond_broadcast(&gs.count_threads_cv);
+        }
+        pthread_mutex_unlock(&gs.count_threads_mutex);
+
+        /* Join exiting threads */
+        for (t=0; t<gs.nthreads; t++) {
+            rc = pthread_join(gs.threads[t], &status);
+            if (rc) {
+                fprintf(stderr,
+                        "ERROR; return code from pthread_join() is %d\n",
+                        rc);
+                fprintf(stderr, "\tError detail: %s\n", strerror(rc));
+                exit(-1);
+            }
+        }
+        gs.init_threads_done = 0;
+        gs.end_threads = 0;
+    }
+
+    /* Launch a new pool of threads (if necessary) */
+    gs.nthreads = nthreads_new;
+    init_threads();
+
+    return nthreads_old;
+}
+
+
+#ifdef USE_VML
+
+static PyObject *
+_get_vml_version(PyObject *self, PyObject *args)
+{
+    int len=198;
+    char buf[198];
+    mkl_get_version_string(buf, len);
+    return Py_BuildValue("s", buf);
+}
+
+static PyObject *
+_set_vml_accuracy_mode(PyObject *self, PyObject *args)
+{
+    int mode_in, mode_old;
+    if (!PyArg_ParseTuple(args, "i", &mode_in))
+    return NULL;
+    mode_old = vmlGetMode() & VML_ACCURACY_MASK;
+    vmlSetMode((mode_in & VML_ACCURACY_MASK) | VML_ERRMODE_IGNORE );
+    return Py_BuildValue("i", mode_old);
+}
+
+static PyObject *
+_set_vml_num_threads(PyObject *self, PyObject *args)
+{
+    int max_num_threads;
+    if (!PyArg_ParseTuple(args, "i", &max_num_threads))
+    return NULL;
+    mkl_domain_set_num_threads(max_num_threads, MKL_DOMAIN_VML);
+    Py_RETURN_NONE;
+}
+
+static PyObject *
+_get_vml_num_threads(PyObject *self, PyObject *args)
+{
+    int max_num_threads = mkl_domain_get_max_threads (MKL_DOMAIN_VML);
+    return Py_BuildValue("i", max_num_threads);
+}
+
+#endif
+
+static PyObject*
+Py_set_num_threads(PyObject *self, PyObject *args)
+{
+    int num_threads, nthreads_old;
+    if (!PyArg_ParseTuple(args, "i", &num_threads))
+    return NULL;
+    nthreads_old = numexpr_set_nthreads(num_threads);
+    return Py_BuildValue("i", nthreads_old);
+}
+
+static PyObject*
+Py_get_num_threads(PyObject *self, PyObject *args)
+{
+    int n_thread;
+    n_thread = gs.nthreads;
+    return Py_BuildValue("i", n_thread);
+}
+
+static PyMethodDef module_methods[] = {
+#ifdef USE_VML
+    {"_get_vml_version", _get_vml_version, METH_VARARGS,
+     "Get the VML/MKL library version."},
+    {"_set_vml_accuracy_mode", _set_vml_accuracy_mode, METH_VARARGS,
+     "Set accuracy mode for VML functions."},
+    {"_set_vml_num_threads", _set_vml_num_threads, METH_VARARGS,
+     "Suggests a maximum number of threads to be used in VML operations."},
+    {"_get_vml_num_threads", _get_vml_num_threads, METH_VARARGS,
+     "Gets the maximum number of threads to be used in VML operations."},
+#endif
+    {"_set_num_threads", Py_set_num_threads, METH_VARARGS,
+     "Suggests a maximum number of threads to be used in operations."},
+    {"_get_num_threads", Py_get_num_threads, METH_VARARGS,
+     "Gets the maximum number of threads currently in use for operations."},
+    {NULL}
+};
+
+static int
+add_symbol(PyObject *d, const char *sname, int name, const char* routine_name)
+{
+    PyObject *o, *s;
+    int r;
+
+    if (!sname) {
+        return 0;
+    }
+
+    o = PyLong_FromLong(name);
+    s = PyBytes_FromString(sname);
+    if (!o || !s) {
+        PyErr_SetString(PyExc_RuntimeError, routine_name);
+        r = -1;
+    }
+    else {
+        r = PyDict_SetItem(d, s, o);
+    }
+    Py_XDECREF(o);
+    Py_XDECREF(s);
+    return r;
+}
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* XXX: handle the "global_state" state via moduledef */
+static struct PyModuleDef moduledef = {
+        PyModuleDef_HEAD_INIT,
+        "interpreter",
+        NULL,
+        -1,                 /* sizeof(struct global_state), */
+        module_methods,
+        NULL,
+        NULL,               /* module_traverse, */
+        NULL,               /* module_clear, */
+        NULL
+};
+
+#define INITERROR return NULL
+
+PyObject *
+PyInit_interpreter(void) {
+    PyObject *m, *d;
+
+
+    char *max_thread_str = getenv("NUMEXPR_MAX_THREADS");
+    char *end;
+    if (max_thread_str != NULL) {
+        global_max_threads = strtol(max_thread_str, &end, 10);
+    }
+
+
+
+    th_params.memsteps    = (npy_intp**)calloc(sizeof(npy_intp*), global_max_threads);
+    th_params.iter        = (NpyIter**)calloc(sizeof(NpyIter*),   global_max_threads);
+    th_params.reduce_iter = (NpyIter**)calloc(sizeof(NpyIter*),   global_max_threads);
+    gs.threads            = (pthread_t*)calloc(sizeof(pthread_t), global_max_threads);
+    gs.tids               = (int*)calloc(sizeof(int),             global_max_threads);
+    // TODO: for Py3, deallocate: https://docs.python.org/3/c-api/module.html#c.PyModuleDef.m_free
+    // For Python 2.7, people have to exit the process to reclaim the memory.
+
+    if (PyType_Ready(&NumExprType) < 0)
+        INITERROR;
+
+    m = PyModule_Create(&moduledef);
+
+    if (m == NULL)
+        INITERROR;
+
+    #ifdef Py_GIL_DISABLED
+        PyUnstable_Module_SetGIL(m, Py_MOD_GIL_NOT_USED);
+    #endif
+
+    Py_INCREF(&NumExprType);
+    PyModule_AddObject(m, "NumExpr", (PyObject *)&NumExprType);
+
+    import_array();
+
+    d = PyDict_New();
+    if (!d) INITERROR;
+
+#define OPCODE(n, name, sname, ...)                              \
+    if (add_symbol(d, sname, name, "add_op") < 0) { INITERROR; }
+#include "opcodes.hpp"
+#undef OPCODE
+
+    if (PyModule_AddObject(m, "opcodes", d) < 0) INITERROR;
+
+    d = PyDict_New();
+    if (!d) INITERROR;
+
+#define add_func(name, sname)                           \
+    if (add_symbol(d, sname, name, "add_func") < 0) { INITERROR; }
+#define FUNC_FF(name, sname, ...)  add_func(name, sname);
+#define FUNC_FFF(name, sname, ...) add_func(name, sname);
+#define FUNC_DD(name, sname, ...)  add_func(name, sname);
+#define FUNC_BF(name, sname, ...)  add_func(name, sname);
+#define FUNC_BD(name, sname, ...)  add_func(name, sname);
+#define FUNC_BC(name, sname, ...)  add_func(name, sname);
+#define FUNC_DDD(name, sname, ...) add_func(name, sname);
+#define FUNC_CC(name, sname, ...)  add_func(name, sname);
+#define FUNC_CCC(name, sname, ...) add_func(name, sname);
+#define FUNC_II(name, sname, ...) add_func(name, sname);
+#define FUNC_LL(name, sname, ...) add_func(name, sname);
+#include "functions.hpp"
+#undef FUNC_LL
+#undef FUNC_II
+#undef FUNC_CCC
+#undef FUNC_CC
+#undef FUNC_DDD
+#undef FUNC_BC
+#undef FUNC_BD
+#undef FUNC_BF
+#undef FUNC_DD
+#undef FUNC_FFF
+#undef FUNC_FF
+#undef add_func
+
+    if (PyModule_AddObject(m, "funccodes", d) < 0) INITERROR;
+
+    if (PyModule_AddObject(m, "allaxes", PyLong_FromLong(255)) < 0) INITERROR;
+    if (PyModule_AddObject(m, "maxdims", PyLong_FromLong(NPY_MAXDIMS)) < 0) INITERROR;
+
+    if(PyModule_AddIntConstant(m, "MAX_THREADS", global_max_threads) < 0) INITERROR;
+
+    // Let's export the block sizes to Python side for benchmarking comparisons
+    if(PyModule_AddIntConstant(m, "__BLOCK_SIZE1__", BLOCK_SIZE1) < 0) INITERROR;
+    // Export if we are using VML or not
+#ifdef USE_VML
+    if(PyModule_AddObject(m, "use_vml", Py_True) < 0) INITERROR;
+#else
+    if(PyModule_AddObject(m, "use_vml", Py_False) < 0) INITERROR;
+#endif
+
+    return m;
+}
+
+#ifdef __cplusplus
+}  // extern "C"
+#endif

numexpr/module.hpp

@@ -0,0 +1,60 @@
+#ifndef NUMEXPR_MODULE_HPP
+#define NUMEXPR_MODULE_HPP
+
+// Deal with the clunky numpy import mechanism
+// by inverting the logic of the NO_IMPORT_ARRAY symbol.
+#define PY_ARRAY_UNIQUE_SYMBOL numexpr_ARRAY_API
+#ifndef DO_NUMPY_IMPORT_ARRAY
+#  define NO_IMPORT_ARRAY
+#endif
+
+#define NPY_NO_DEPRECATED_API NPY_API_VERSION
+
+#include <Python.h>
+#include <numpy/ndarrayobject.h>
+#include <numpy/arrayscalars.h>
+
+#include "numexpr_config.hpp"
+
+struct global_state {
+    /* Global variables for threads */
+    int nthreads;                    /* number of desired threads in pool */
+    int init_threads_done;           /* pool of threads initialized? */
+    int end_threads;                 /* should exisiting threads end? */
+    // pthread_t threads[MAX_THREADS];  /* opaque structure for threads */
+    // int tids[MAX_THREADS];           /* ID per each thread */
+    /* NOTE: threads and tids are arrays, they MUST be allocated to length
+       `global_max_threads` before module load. */
+    pthread_t *threads;              /* opaque structure for threads */
+    int *tids;                       /* ID per each thread */
+    npy_intp gindex;                 /* global index for all threads */
+    int init_sentinels_done;         /* sentinels initialized? */
+    int giveup;                      /* should parallel code giveup? */
+    int force_serial;                /* force serial code instead of parallel? */
+    int pid;                         /* the PID for this process */
+
+    /* Synchronization variables for threadpool state */
+    pthread_mutex_t count_mutex;
+    int count_threads;
+    int barrier_passed;         /* indicates if the thread pool's thread barrier
+                                   is unlocked and ready for the VM to process.*/
+    pthread_mutex_t count_threads_mutex;
+    pthread_cond_t count_threads_cv;
+
+    /* Mutual exclusion for access to global thread params (th_params) */
+    pthread_mutex_t parallel_mutex;
+
+    global_state() {
+        nthreads = 1;
+        init_threads_done = 0;
+        barrier_passed = 0;
+        end_threads = 0;
+        pid = 0;
+    }
+};
+
+extern global_state gs;
+
+int numexpr_set_nthreads(int nthreads_new);
+
+#endif // NUMEXPR_MODULE_HPP

numexpr/msvc_function_stubs.hpp

@@ -0,0 +1,231 @@
+#include <float.h>  // for _finite, _isnan on MSVC
+
+#ifndef NUMEXPR_MSVC_FUNCTION_STUBS_HPP
+#define NUMEXPR_MSVC_FUNCTION_STUBS_HPP
+
+/*********************************************************************
+  Numexpr - Fast numerical array expression evaluator for NumPy.
+
+      License: MIT
+      Author:  See AUTHORS.txt
+
+  See LICENSE.txt for details about copyright and rights to use.
+**********************************************************************/
+
+/* Declare stub functions for MSVC.  It turns out that single precision
+   definitions in <math.h> are actually #define'd and are not usable
+   as function pointers :-/ */
+
+/* Due to casting problems (normally return ints not bools, easiest to define
+non-overloaded wrappers for these functions) */
+// MSVC version: use global ::isfinite / ::isnan
+inline bool isfinitef_(float x) { return !!::_finite(x); }   // MSVC has _finite
+inline bool isnanf_(float x)    { return !!::_isnan(x); }    // MSVC has _isnan
+inline bool isfinited(double x) { return !!::_finite(x); }
+inline bool isnand(double x)    { return !!::_isnan(x); }
+inline bool isinfd(double x) { return !!::isinf(x); }
+inline bool isinff_(float x)    { return !!::isinf(x); }
+
+// To handle overloading of fmax/fmin in cmath and match NumPy behaviour for NaNs
+inline double fmaxd(double x, double y)    { return (isnand(x) | isnand(y))? NAN : fmax(x, y); }
+inline double fmind(double x, double y)    { return (isnand(x) | isnand(y))? NAN : fmin(x, y); }
+
+
+#if _MSC_VER < 1400  // 1310 == MSVC 7.1
+    /* Apparently, single precision functions are not included in MSVC 7.1 */
+
+    #define sqrtf(x)    ((float)sqrt((double)(x)))
+    #define sinf(x)    ((float)sin((double)(x)))
+    #define cosf(x)    ((float)cos((double)(x)))
+    #define tanf(x)    ((float)tan((double)(x)))
+    #define asinf(x)    ((float)asin((double)(x)))
+    #define acosf(x)    ((float)acos((double)(x)))
+    #define atanf(x)    ((float)atan((double)(x)))
+    #define sinhf(x)    ((float)sinh((double)(x)))
+    #define coshf(x)    ((float)cosh((double)(x)))
+    #define tanhf(x)    ((float)tanh((double)(x)))
+    #define asinhf(x)    ((float)asinh((double)(x)))
+    #define acoshf(x)    ((float)acosh((double)(x)))
+    #define atanhf(x)    ((float)atanh((double)(x)))
+    #define logf(x)    ((float)log((double)(x)))
+    #define log1pf(x)    ((float)log1p((double)(x)))
+    #define log10f(x)    ((float)log10((double)(x)))
+    #define log2f(x)    ((float)log2((double)(x)))
+    #define expf(x)    ((float)exp((double)(x)))
+    #define expm1f(x)    ((float)expm1((double)(x)))
+    #define fabsf(x)    ((float)fabs((double)(x)))
+    #define fmodf(x, y)    ((float)fmod((double)(x), (double)(y)))
+    #define atan2f(x, y)    ((float)atan2((double)(x), (double)(y)))
+    #define hypotf(x, y)    ((float)hypot((double)(x), (double)(y)))
+    #define copysignf(x, y)    ((float)copysign((double)(x), (double)(y)))
+    #define nextafterf(x, y)    ((float)nextafter((double)(x), (double)(y)))
+    #define ceilf(x)    ((float)ceil((double)(x)))
+    #define hypotf(x)    ((float)hypot((double)(x)))
+    #define rintf(x)    ((float)rint((double)(x)))
+    #define truncf(x)    ((float)trunc((double)(x)))
+
+
+    /* The next are directly called from interp_body.cpp */
+    #define powf(x, y)    ((float)pow((double)(x), (double)(y)))
+    #define floorf(x)    ((float)floor((double)(x)))
+
+    #define fmaxf_(x, y)    ((float)fmaxd((double)(x), (double)(y))) // define fmaxf_ since fmaxf doesn't exist for early MSVC
+    #define fminf_(x, y)    ((float)fmind((double)(x), (double)(y)))
+#else
+    inline float fmaxf_(float x, float y)    { return (isnanf_(x) | isnanf_(y))? NAN : fmaxf(x, y); }
+    inline float fminf_(float x, float y)    { return (isnanf_(x) | isnanf_(y))? NAN : fminf(x, y); }
+#endif // _MSC_VER < 1400
+
+
+/* Now the actual stubs */
+
+inline float sqrtf2(float x) {
+    return sqrtf(x);
+}
+
+inline float sinf2(float x) {
+    return sinf(x);
+}
+
+inline float cosf2(float x) {
+    return cosf(x);
+}
+
+inline float tanf2(float x) {
+    return tanf(x);
+}
+
+inline float asinf2(float x) {
+    return asinf(x);
+}
+
+inline float acosf2(float x) {
+    return acosf(x);
+}
+
+inline float atanf2(float x) {
+    return atanf(x);
+}
+
+inline float sinhf2(float x) {
+    return sinhf(x);
+}
+
+inline float coshf2(float x) {
+    return coshf(x);
+}
+
+inline float tanhf2(float x) {
+    return tanhf(x);
+}
+
+inline float asinhf2(float x) {
+    return asinhf(x);
+}
+
+inline float acoshf2(float x) {
+    return acoshf(x);
+}
+
+inline float atanhf2(float x) {
+    return atanhf(x);
+}
+
+inline float logf2(float x) {
+    return logf(x);
+}
+
+inline float log1pf2(float x) {
+    return log1pf(x);
+}
+
+inline float log10f2(float x) {
+    return log10f(x);
+}
+
+inline float log2f2(float x) {
+    return log2f(x);
+}
+
+inline float expf2(float x) {
+    return expf(x);
+}
+
+inline float expm1f2(float x) {
+    return expm1f(x);
+}
+
+inline float fabsf2(float x) {
+    return fabsf(x);
+}
+
+inline float fmodf2(float x, float y) {
+    return fmodf(x, y);
+}
+
+inline float atan2f2(float x, float y) {
+    return atan2f(x, y);
+}
+
+inline float hypotf2(float x, float y) {
+    return hypotf(x, y);
+}
+
+inline float nextafterf2(float x, float y) {
+    return nextafterf(x, y);
+}
+
+inline float copysignf2(float x, float y) {
+    return copysignf(x, y);
+}
+
+inline float fmaxf2(float x, float y) {
+    return fmaxf_(x, y);
+}
+
+inline float fminf2(float x, float y) {
+    return fminf_(x, y);
+}
+
+
+// Boolean output functions
+inline bool isnanf2(float x) {
+    return isnanf_(x);
+}
+
+inline bool isfinitef2(float x) {
+    return isfinitef_(x);
+}
+
+inline bool isinff2(float x) {
+    return isinff_(x);
+}
+
+
+// Needed for allowing the internal casting in numexpr machinery for
+// conjugate operations
+inline float fconjf2(float x) {
+    return x;
+}
+
+inline float ceilf2(float x) {
+    return ceilf(x);
+}
+
+inline float floorf2(float x) {
+    return floorf(x);
+}
+
+inline float rintf2(float x) {
+    return rintf(x);
+}
+
+inline float truncf2(float x) {
+    return truncf(x);
+}
+
+inline bool signbitf2(float x) {
+    return signbitf(x);
+}
+
+#endif // NUMEXPR_MSVC_FUNCTION_STUBS_HPP

numexpr/numexpr_config.hpp

@@ -0,0 +1,71 @@
+#ifndef NUMEXPR_CONFIG_HPP
+#define NUMEXPR_CONFIG_HPP
+
+// x86 platform works with unaligned reads and writes
+// MW: I have seen exceptions to this when the compiler chooses to use aligned SSE
+#if (defined(NPY_CPU_X86) || defined(NPY_CPU_AMD64))
+#  define USE_UNALIGNED_ACCESS 1
+#endif
+
+// #ifdef SCIPY_MKL_H
+// #define USE_VML
+// #endif
+
+#ifdef USE_VML
+/* The values below have been tuned for a Skylake processor (E3-1245 v5 @ 3.50GHz) */
+#define BLOCK_SIZE1 1024
+#else
+/* The values below have been tuned for a Skylake processor (E3-1245 v5 @ 3.50GHz) */
+#define BLOCK_SIZE1 1024
+#endif
+
+// The default threadpool size. It's prefer that the user set this via an
+// environment variable, "NUMEXPR_MAX_THREADS"
+#define DEFAULT_MAX_THREADS 64
+
+// Remove dependence on NPY_MAXARGS, which would be a runtime constant instead of compiletime
+// constant. If numpy raises NPY_MAXARGS, we should notice and raise this as well
+#define NE_MAXARGS 64
+
+#if defined(_WIN32)
+  #include "win32/pthread.h"
+  #include <process.h>
+  #define getpid _getpid
+#else
+  #include <pthread.h>
+  #include "unistd.h"
+#endif
+
+#ifdef USE_VML
+#include "mkl_vml.h"
+#include "mkl_service.h"
+#endif
+#include <cmath>
+//no single precision version of signbit in C++ standard
+inline bool signbitf(float x) { return signbit((double)x); }
+
+#ifdef _WIN32
+  #ifndef __MINGW32__
+    #include "missing_posix_functions.hpp"
+  #endif
+  #include "msvc_function_stubs.hpp"
+#else
+/* GCC/Clang version: use std:: (can't use it for windows)
+   msvc_function_stubs contains windows alternatives */
+/* Due to casting problems (normally return ints not bools, easiest to define
+   non-overloaded wrappers for these functions) */
+inline bool isfinitef_(float x) { return !!std::isfinite(x); }
+inline bool isnanf_(float x)    { return !!std::isnan(x); }
+inline bool isfinited(double x) { return !!std::isfinite(x); }
+inline bool isnand(double x)    { return !!std::isnan(x); }
+inline bool isinff_(float x) { return !!std::isinf(x); }
+inline bool isinfd(double x)    { return !!std::isinf(x); }
+
+// To handle overloading of fmax/fmin in cmath and match NumPy behaviour for NaNs
+inline double fmaxd(double x, double y)    { return (isnand(x) | isnand(y))? NAN : fmax(x, y); }
+inline double fmind(double x, double y)    { return (isnand(x) | isnand(y))? NAN : fmin(x, y); }
+inline float fmaxf_(float x, float y)    { return (isnanf_(x) | isnanf_(y))? NAN : fmaxf(x, y); }
+inline float fminf_(float x, float y)    { return (isnanf_(x) | isnanf_(y))? NAN : fminf(x, y); }
+#endif
+
+#endif // NUMEXPR_CONFIG_HPP

numexpr/numexpr_object.cpp

@@ -0,0 +1,407 @@
+/*********************************************************************
+  Numexpr - Fast numerical array expression evaluator for NumPy.
+
+      License: MIT
+      Author:  See AUTHORS.txt
+
+  See LICENSE.txt for details about copyright and rights to use.
+**********************************************************************/
+
+#include "module.hpp"
+#include <structmember.h>
+
+#include "numexpr_config.hpp"
+#include "interpreter.hpp"
+#include "numexpr_object.hpp"
+
+static int
+size_from_char(char c)
+{
+    switch (c) {
+        case 'b': return sizeof(char);
+        case 'i': return sizeof(int);
+        case 'l': return sizeof(long long);
+        case 'f': return sizeof(float);
+        case 'd': return sizeof(double);
+        case 'c': return 2*sizeof(double);
+        case 's': return 0;  /* strings are ok but size must be computed */
+        default:
+            PyErr_SetString(PyExc_TypeError, "signature value not in 'bilfdcs'");
+            return -1;
+    }
+}
+
+static void
+NumExpr_dealloc(NumExprObject *self)
+{
+    Py_XDECREF(self->signature);
+    Py_XDECREF(self->tempsig);
+    Py_XDECREF(self->constsig);
+    Py_XDECREF(self->fullsig);
+    Py_XDECREF(self->program);
+    Py_XDECREF(self->constants);
+    Py_XDECREF(self->input_names);
+    PyMem_Del(self->mem);
+    PyMem_Del(self->rawmem);
+    PyMem_Del(self->memsteps);
+    PyMem_Del(self->memsizes);
+    Py_TYPE(self)->tp_free((PyObject*)self);
+}
+
+static PyObject *
+NumExpr_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    NumExprObject *self = (NumExprObject *)type->tp_alloc(type, 0);
+    if (self != NULL) {
+#define INIT_WITH(name, object) \
+        self->name = object; \
+        if (!self->name) { \
+            Py_DECREF(self); \
+            return NULL; \
+        }
+
+        INIT_WITH(signature, PyBytes_FromString(""));
+        INIT_WITH(tempsig, PyBytes_FromString(""));
+        INIT_WITH(constsig, PyBytes_FromString(""));
+        INIT_WITH(fullsig, PyBytes_FromString(""));
+        INIT_WITH(program, PyBytes_FromString(""));
+        INIT_WITH(constants, PyTuple_New(0));
+        Py_INCREF(Py_None);
+        self->input_names = Py_None;
+        self->mem = NULL;
+        self->rawmem = NULL;
+        self->memsteps = NULL;
+        self->memsizes = NULL;
+        self->rawmemsize = 0;
+        self->n_inputs = 0;
+        self->n_constants = 0;
+        self->n_temps = 0;
+#undef INIT_WITH
+    }
+    return (PyObject *)self;
+}
+
+#define CHARP(s) ((char *)(s))
+
+static int
+NumExpr_init(NumExprObject *self, PyObject *args, PyObject *kwds)
+{
+    int i, j, mem_offset;
+    int n_inputs, n_constants, n_temps;
+    PyObject *signature = NULL, *tempsig = NULL, *constsig = NULL;
+    PyObject *fullsig = NULL, *program = NULL, *constants = NULL;
+    PyObject *input_names = NULL, *o_constants = NULL;
+    int *itemsizes = NULL;
+    char **mem = NULL, *rawmem = NULL;
+    npy_intp *memsteps;
+    npy_intp *memsizes;
+    int rawmemsize;
+    static char *kwlist[] = {CHARP("signature"), CHARP("tempsig"),
+			     CHARP("program"),  CHARP("constants"),
+			     CHARP("input_names"), NULL};
+
+    if (!PyArg_ParseTupleAndKeywords(args, kwds, "SSS|OO", kwlist,
+                                     &signature,
+                                     &tempsig,
+                                     &program, &o_constants,
+                                     &input_names)) {
+        return -1;
+    }
+
+    n_inputs = (int)PyBytes_Size(signature);
+    n_temps = (int)PyBytes_Size(tempsig);
+
+    if (o_constants) {
+        if (!PySequence_Check(o_constants) ) {
+                PyErr_SetString(PyExc_TypeError, "constants must be a sequence");
+                return -1;
+        }
+        n_constants = (int)PySequence_Length(o_constants);
+        if (!(constants = PyTuple_New(n_constants)))
+            return -1;
+        if (!(constsig = PyBytes_FromStringAndSize(NULL, n_constants))) {
+            Py_DECREF(constants);
+            return -1;
+        }
+        if (!(itemsizes = PyMem_New(int, n_constants))) {
+            Py_DECREF(constants);
+            return -1;
+        }
+        for (i = 0; i < n_constants; i++) {
+            PyObject *o;
+            if (!(o = PySequence_GetItem(o_constants, i))) { /* new reference */
+                Py_DECREF(constants);
+                Py_DECREF(constsig);
+                PyMem_Del(itemsizes);
+                return -1;
+            }
+            PyTuple_SET_ITEM(constants, i, o); /* steals reference */
+            if (PyBool_Check(o)) {
+                PyBytes_AS_STRING(constsig)[i] = 'b';
+                itemsizes[i] = size_from_char('b');
+                continue;
+            }
+
+            if (PyArray_IsScalar(o, Int32)) {
+                PyBytes_AS_STRING(constsig)[i] = 'i';
+                itemsizes[i] = size_from_char('i');
+                continue;
+            }
+
+            if (PyArray_IsScalar(o, Int64)) {
+                PyBytes_AS_STRING(constsig)[i] = 'l';
+                itemsizes[i] = size_from_char('l');
+                continue;
+            }
+            /* The Float32 scalars are the only ones that should reach here */
+            if (PyArray_IsScalar(o, Float32)) {
+                PyBytes_AS_STRING(constsig)[i] = 'f';
+                itemsizes[i] = size_from_char('f');
+                continue;
+            }
+            if (PyFloat_Check(o)) {
+                /* Python float constants are double precision by default */
+                PyBytes_AS_STRING(constsig)[i] = 'd';
+                itemsizes[i] = size_from_char('d');
+                continue;
+            }
+            if (PyComplex_Check(o)) {
+                PyBytes_AS_STRING(constsig)[i] = 'c';
+                itemsizes[i] = size_from_char('c');
+                continue;
+            }
+            if (PyBytes_Check(o)) {
+                PyBytes_AS_STRING(constsig)[i] = 's';
+                itemsizes[i] = (int)PyBytes_GET_SIZE(o);
+                continue;
+            }
+            PyErr_SetString(PyExc_TypeError, "constants must be of type bool/int/long/float/double/complex/bytes");
+            Py_DECREF(constsig);
+            Py_DECREF(constants);
+            PyMem_Del(itemsizes);
+            return -1;
+        }
+    } else {
+        n_constants = 0;
+        if (!(constants = PyTuple_New(0)))
+            return -1;
+        if (!(constsig = PyBytes_FromString(""))) {
+            Py_DECREF(constants);
+            return -1;
+        }
+    }
+
+    fullsig = PyBytes_FromFormat("%c%s%s%s", get_return_sig(program),
+        PyBytes_AS_STRING(signature), PyBytes_AS_STRING(constsig),
+        PyBytes_AS_STRING(tempsig));
+    if (!fullsig) {
+        Py_DECREF(constants);
+        Py_DECREF(constsig);
+        PyMem_Del(itemsizes);
+        return -1;
+    }
+
+    if (!input_names) {
+        input_names = Py_None;
+    }
+
+    /* Compute the size of registers. We leave temps out (will be
+       malloc'ed later on). */
+    rawmemsize = 0;
+    for (i = 0; i < n_constants; i++)
+        rawmemsize += itemsizes[i];
+    rawmemsize *= BLOCK_SIZE1;
+
+    mem = PyMem_New(char *, 1 + n_inputs + n_constants + n_temps);
+    rawmem = PyMem_New(char, rawmemsize);
+    memsteps = PyMem_New(npy_intp, 1 + n_inputs + n_constants + n_temps);
+    memsizes = PyMem_New(npy_intp, 1 + n_inputs + n_constants + n_temps);
+    if (!mem || !rawmem || !memsteps || !memsizes) {
+        Py_DECREF(constants);
+        Py_DECREF(constsig);
+        Py_DECREF(fullsig);
+        PyMem_Del(itemsizes);
+        PyMem_Del(mem);
+        PyMem_Del(rawmem);
+        PyMem_Del(memsteps);
+        PyMem_Del(memsizes);
+        return -1;
+    }
+    /*
+       0                                                  -> output
+       [1, n_inputs+1)                                    -> inputs
+       [n_inputs+1, n_inputs+n_consts+1)                  -> constants
+       [n_inputs+n_consts+1, n_inputs+n_consts+n_temps+1) -> temps
+    */
+    /* Fill in 'mem' and 'rawmem' for constants */
+    mem_offset = 0;
+    for (i = 0; i < n_constants; i++) {
+        char c = PyBytes_AS_STRING(constsig)[i];
+        int size = itemsizes[i];
+        mem[i+n_inputs+1] = rawmem + mem_offset;
+        mem_offset += BLOCK_SIZE1 * size;
+        memsteps[i+n_inputs+1] = memsizes[i+n_inputs+1] = size;
+        /* fill in the constants */
+        if (c == 'b') {
+            char *bmem = (char*)mem[i+n_inputs+1];
+            char value = (char)PyLong_AsLong(PyTuple_GET_ITEM(constants, i));
+            for (j = 0; j < BLOCK_SIZE1; j++) {
+                bmem[j] = value;
+            }
+        } else if (c == 'i') {
+            int *imem = (int*)mem[i+n_inputs+1];
+            int value = (int)PyLong_AsLong(PyTuple_GET_ITEM(constants, i));
+            for (j = 0; j < BLOCK_SIZE1; j++) {
+                imem[j] = value;
+            }
+        } else if (c == 'l') {
+            long long *lmem = (long long*)mem[i+n_inputs+1];
+            long long value = PyLong_AsLongLong(PyTuple_GET_ITEM(constants, i));
+            for (j = 0; j < BLOCK_SIZE1; j++) {
+                lmem[j] = value;
+            }
+        } else if (c == 'f') {
+            /* In this particular case the constant is in a NumPy scalar
+             and in a regular Python object */
+            float *fmem = (float*)mem[i+n_inputs+1];
+            float value = PyArrayScalar_VAL(PyTuple_GET_ITEM(constants, i),
+                                            Float);
+            for (j = 0; j < BLOCK_SIZE1; j++) {
+                fmem[j] = value;
+            }
+        } else if (c == 'd') {
+            double *dmem = (double*)mem[i+n_inputs+1];
+            double value = PyFloat_AS_DOUBLE(PyTuple_GET_ITEM(constants, i));
+            for (j = 0; j < BLOCK_SIZE1; j++) {
+                dmem[j] = value;
+            }
+        } else if (c == 'c') {
+            double *cmem = (double*)mem[i+n_inputs+1];
+            Py_complex value = PyComplex_AsCComplex(PyTuple_GET_ITEM(constants, i));
+            for (j = 0; j < 2*BLOCK_SIZE1; j+=2) {
+                cmem[j] = value.real;
+                cmem[j+1] = value.imag;
+            }
+        } else if (c == 's') {
+            char *smem = (char*)mem[i+n_inputs+1];
+            char *value = PyBytes_AS_STRING(PyTuple_GET_ITEM(constants, i));
+            for (j = 0; j < size*BLOCK_SIZE1; j+=size) {
+                memcpy(smem + j, value, size);
+            }
+        }
+    }
+    /* This is no longer needed since no unusual item sizes appear
+       in temporaries (there are no string temporaries). */
+    PyMem_Del(itemsizes);
+
+    /* Fill in 'memsteps' and 'memsizes' for temps */
+    for (i = 0; i < n_temps; i++) {
+        char c = PyBytes_AS_STRING(tempsig)[i];
+        int size = size_from_char(c);
+        memsteps[i+n_inputs+n_constants+1] = size;
+        memsizes[i+n_inputs+n_constants+1] = size;
+    }
+    /* See if any errors occured (e.g., in size_from_char) or if mem_offset is wrong */
+    if (PyErr_Occurred() || mem_offset != rawmemsize) {
+        if (mem_offset != rawmemsize) {
+            PyErr_Format(PyExc_RuntimeError, "mem_offset does not match rawmemsize");
+        }
+        Py_DECREF(constants);
+        Py_DECREF(constsig);
+        Py_DECREF(fullsig);
+        PyMem_Del(mem);
+        PyMem_Del(rawmem);
+        PyMem_Del(memsteps);
+        PyMem_Del(memsizes);
+        return -1;
+    }
+
+
+    #define REPLACE_OBJ(arg) \
+    {PyObject *tmp = self->arg; \
+     self->arg = arg; \
+     Py_XDECREF(tmp);}
+    #define INCREF_REPLACE_OBJ(arg) {Py_INCREF(arg); REPLACE_OBJ(arg);}
+    #define REPLACE_MEM(arg) {PyMem_Del(self->arg); self->arg=arg;}
+
+    INCREF_REPLACE_OBJ(signature);
+    INCREF_REPLACE_OBJ(tempsig);
+    REPLACE_OBJ(constsig);
+    REPLACE_OBJ(fullsig);
+    INCREF_REPLACE_OBJ(program);
+    REPLACE_OBJ(constants);
+    INCREF_REPLACE_OBJ(input_names);
+    REPLACE_MEM(mem);
+    REPLACE_MEM(rawmem);
+    REPLACE_MEM(memsteps);
+    REPLACE_MEM(memsizes);
+    self->rawmemsize = rawmemsize;
+    self->n_inputs = n_inputs;
+    self->n_constants = n_constants;
+    self->n_temps = n_temps;
+
+    #undef REPLACE_OBJ
+    #undef INCREF_REPLACE_OBJ
+    #undef REPLACE_MEM
+
+    return check_program(self);
+}
+
+static PyMethodDef NumExpr_methods[] = {
+    {"run", (PyCFunction) NumExpr_run, METH_VARARGS|METH_KEYWORDS, NULL},
+    {NULL, NULL}
+};
+
+static PyMemberDef NumExpr_members[] = {
+    {CHARP("signature"), T_OBJECT_EX, offsetof(NumExprObject, signature), READONLY, NULL},
+    {CHARP("constsig"), T_OBJECT_EX, offsetof(NumExprObject, constsig), READONLY, NULL},
+    {CHARP("tempsig"), T_OBJECT_EX, offsetof(NumExprObject, tempsig), READONLY, NULL},
+    {CHARP("fullsig"), T_OBJECT_EX, offsetof(NumExprObject, fullsig), READONLY, NULL},
+
+    {CHARP("program"), T_OBJECT_EX, offsetof(NumExprObject, program), READONLY, NULL},
+    {CHARP("constants"), T_OBJECT_EX, offsetof(NumExprObject, constants),
+     READONLY, NULL},
+    {CHARP("input_names"), T_OBJECT, offsetof(NumExprObject, input_names), 0, NULL},
+    {NULL},
+};
+
+PyTypeObject NumExprType = {
+    PyVarObject_HEAD_INIT(NULL, 0)
+    "numexpr.NumExpr",         /*tp_name*/
+    sizeof(NumExprObject),     /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)NumExpr_dealloc, /*tp_dealloc*/
+    0,                         /*tp_print*/
+    0,                         /*tp_getattr*/
+    0,                         /*tp_setattr*/
+    0,                         /*tp_compare*/
+    0,                         /*tp_repr*/
+    0,                         /*tp_as_number*/
+    0,                         /*tp_as_sequence*/
+    0,                         /*tp_as_mapping*/
+    0,                         /*tp_hash */
+    (ternaryfunc)NumExpr_run,  /*tp_call*/
+    0,                         /*tp_str*/
+    0,                         /*tp_getattro*/
+    0,                         /*tp_setattro*/
+    0,                         /*tp_as_buffer*/
+    Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
+    "NumExpr objects",         /* tp_doc */
+    0,                       /* tp_traverse */
+    0,                       /* tp_clear */
+    0,                       /* tp_richcompare */
+    0,                       /* tp_weaklistoffset */
+    0,                       /* tp_iter */
+    0,                       /* tp_iternext */
+    NumExpr_methods,           /* tp_methods */
+    NumExpr_members,           /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)NumExpr_init,    /* tp_init */
+    0,                         /* tp_alloc */
+    NumExpr_new,               /* tp_new */
+};

numexpr/numexpr_object.hpp

@@ -0,0 +1,34 @@
+#ifndef NUMEXPR_OBJECT_HPP
+#define NUMEXPR_OBJECT_HPP
+/*********************************************************************
+  Numexpr - Fast numerical array expression evaluator for NumPy.
+
+      License: MIT
+      Author:  See AUTHORS.txt
+
+  See LICENSE.txt for details about copyright and rights to use.
+**********************************************************************/
+
+struct NumExprObject
+{
+    PyObject_HEAD
+    PyObject *signature;    /* a python string */
+    PyObject *tempsig;
+    PyObject *constsig;
+    PyObject *fullsig;
+    PyObject *program;      /* a python string */
+    PyObject *constants;    /* a tuple of int/float/complex */
+    PyObject *input_names;  /* tuple of strings */
+    char **mem;             /* pointers to registers */
+    char *rawmem;           /* a chunks of raw memory for storing registers */
+    npy_intp *memsteps;
+    npy_intp *memsizes;
+    int  rawmemsize;
+    int  n_inputs;
+    int  n_constants;
+    int  n_temps;
+};
+
+extern PyTypeObject NumExprType;
+
+#endif // NUMEXPR_OBJECT_HPP

numexpr/opcodes.hpp

@@ -0,0 +1,214 @@
+/*********************************************************************
+  Numexpr - Fast numerical array expression evaluator for NumPy.
+
+      License: MIT
+      Author:  See AUTHORS.txt
+
+  See LICENSE.txt for details about copyright and rights to use.
+**********************************************************************/
+
+/*
+OPCODE(n, enum_name, exported, return_type, arg1_type, arg2_type, arg3_type)
+
+`exported` is NULL if the opcode shouldn't exported by the Python module.
+
+Types are Tb, Ti, Tl, Tf, Td, Tc, Ts, Tn, and T0; these symbols should be
+#defined to whatever is needed. (T0 is the no-such-arg type.)
+
+When adding new OPCODES, one has to respect the order of the numeration, as
+there are parts of the code (iterations) which assume that the OPCODES are ordered.
+
+*/
+OPCODE(0, OP_NOOP, "noop", T0, T0, T0, T0)
+
+OPCODE(1, OP_COPY_BB, "copy_bb", Tb, Tb, T0, T0)
+
+OPCODE(2, OP_INVERT_BB, "invert_bb", Tb, Tb, T0, T0)
+OPCODE(3, OP_AND_BBB, "and_bbb", Tb, Tb, Tb, T0)
+OPCODE(4, OP_OR_BBB, "or_bbb", Tb, Tb, Tb, T0)
+OPCODE(5, OP_XOR_BBB, "xor_bbb", Tb, Tb, Tb, T0)
+
+OPCODE(6, OP_EQ_BBB, "eq_bbb", Tb, Tb, Tb, T0)
+OPCODE(7, OP_NE_BBB, "ne_bbb", Tb, Tb, Tb, T0)
+
+OPCODE(8, OP_GT_BII, "gt_bii", Tb, Ti, Ti, T0)
+OPCODE(9, OP_GE_BII, "ge_bii", Tb, Ti, Ti, T0)
+OPCODE(10, OP_EQ_BII, "eq_bii", Tb, Ti, Ti, T0)
+OPCODE(11, OP_NE_BII, "ne_bii", Tb, Ti, Ti, T0)
+
+OPCODE(12, OP_GT_BLL, "gt_bll", Tb, Tl, Tl, T0)
+OPCODE(13, OP_GE_BLL, "ge_bll", Tb, Tl, Tl, T0)
+OPCODE(14, OP_EQ_BLL, "eq_bll", Tb, Tl, Tl, T0)
+OPCODE(15, OP_NE_BLL, "ne_bll", Tb, Tl, Tl, T0)
+
+OPCODE(16, OP_GT_BFF, "gt_bff", Tb, Tf, Tf, T0)
+OPCODE(17, OP_GE_BFF, "ge_bff", Tb, Tf, Tf, T0)
+OPCODE(18, OP_EQ_BFF, "eq_bff", Tb, Tf, Tf, T0)
+OPCODE(19, OP_NE_BFF, "ne_bff", Tb, Tf, Tf, T0)
+
+OPCODE(20, OP_GT_BDD, "gt_bdd", Tb, Td, Td, T0)
+OPCODE(21, OP_GE_BDD, "ge_bdd", Tb, Td, Td, T0)
+OPCODE(22, OP_EQ_BDD, "eq_bdd", Tb, Td, Td, T0)
+OPCODE(23, OP_NE_BDD, "ne_bdd", Tb, Td, Td, T0)
+
+OPCODE(24, OP_GT_BSS, "gt_bss", Tb, Ts, Ts, T0)
+OPCODE(25, OP_GE_BSS, "ge_bss", Tb, Ts, Ts, T0)
+OPCODE(26, OP_EQ_BSS, "eq_bss", Tb, Ts, Ts, T0)
+OPCODE(27, OP_NE_BSS, "ne_bss", Tb, Ts, Ts, T0)
+
+OPCODE(28, OP_CAST_IB, "cast_ib", Ti, Tb, T0, T0)
+OPCODE(29, OP_COPY_II, "copy_ii", Ti, Ti, T0, T0)
+OPCODE(30, OP_ONES_LIKE_II, "ones_like_ii", Ti, T0, T0, T0)
+OPCODE(31, OP_NEG_II, "neg_ii", Ti, Ti, T0, T0)
+OPCODE(32, OP_ADD_III, "add_iii", Ti, Ti, Ti, T0)
+OPCODE(33, OP_SUB_III, "sub_iii", Ti, Ti, Ti, T0)
+OPCODE(34, OP_MUL_III, "mul_iii", Ti, Ti, Ti, T0)
+OPCODE(35, OP_DIV_III, "div_iii", Ti, Ti, Ti, T0)
+OPCODE(36, OP_POW_III, "pow_iii", Ti, Ti, Ti, T0)
+OPCODE(37, OP_MOD_III, "mod_iii", Ti, Ti, Ti, T0)
+OPCODE(38, OP_FLOORDIV_III, "floordiv_iii", Ti, Ti, Ti, T0)
+
+
+OPCODE(39, OP_LSHIFT_III, "lshift_iii", Ti, Ti, Ti, T0)
+OPCODE(40, OP_RSHIFT_III, "rshift_iii", Ti, Ti, Ti, T0)
+
+OPCODE(41, OP_WHERE_IBII, "where_ibii", Ti, Tb, Ti, Ti)
+// Bitwise ops
+OPCODE(42, OP_INVERT_II, "invert_ii", Ti, Ti, T0, T0)
+OPCODE(43, OP_AND_III, "and_iii", Ti, Ti, Ti, T0)
+OPCODE(44, OP_OR_III, "or_iii", Ti, Ti, Ti, T0)
+OPCODE(45, OP_XOR_III, "xor_iii", Ti, Ti, Ti, T0)
+
+OPCODE(46, OP_CAST_LI, "cast_li", Tl, Ti, T0, T0)
+OPCODE(47, OP_COPY_LL, "copy_ll", Tl, Tl, T0, T0)
+OPCODE(48, OP_ONES_LIKE_LL, "ones_like_ll", Tl, T0, T0, T0)
+OPCODE(49, OP_NEG_LL, "neg_ll", Tl, Tl, T0, T0)
+OPCODE(50, OP_ADD_LLL, "add_lll", Tl, Tl, Tl, T0)
+OPCODE(51, OP_SUB_LLL, "sub_lll", Tl, Tl, Tl, T0)
+OPCODE(52, OP_MUL_LLL, "mul_lll", Tl, Tl, Tl, T0)
+OPCODE(53, OP_DIV_LLL, "div_lll", Tl, Tl, Tl, T0)
+OPCODE(54, OP_POW_LLL, "pow_lll", Tl, Tl, Tl, T0)
+OPCODE(55, OP_MOD_LLL, "mod_lll", Tl, Tl, Tl, T0)
+OPCODE(56, OP_FLOORDIV_LLL, "floordiv_lll", Tl, Tl, Tl, T0)
+
+OPCODE(57, OP_LSHIFT_LLL, "lshift_lll", Tl, Tl, Tl, T0)
+OPCODE(58, OP_RSHIFT_LLL, "rshift_lll", Tl, Tl, Tl, T0)
+
+OPCODE(59, OP_WHERE_LBLL, "where_lbll", Tl, Tb, Tl, Tl)
+// Bitwise ops
+OPCODE(60, OP_INVERT_LL, "invert_ll", Tl, Tl, T0, T0)
+OPCODE(61, OP_AND_LLL, "and_lll", Tl, Tl, Tl, T0)
+OPCODE(62, OP_OR_LLL, "or_lll", Tl, Tl, Tl, T0)
+OPCODE(63, OP_XOR_LLL, "xor_lll", Tl, Tl, Tl, T0)
+
+OPCODE(64, OP_CAST_FI, "cast_fi", Tf, Ti, T0, T0)
+OPCODE(65, OP_CAST_FL, "cast_fl", Tf, Tl, T0, T0)
+OPCODE(66, OP_COPY_FF, "copy_ff", Tf, Tf, T0, T0)
+OPCODE(67, OP_ONES_LIKE_FF, "ones_like_ff", Tf, T0, T0, T0)
+OPCODE(68, OP_NEG_FF, "neg_ff", Tf, Tf, T0, T0)
+OPCODE(69, OP_ADD_FFF, "add_fff", Tf, Tf, Tf, T0)
+OPCODE(70, OP_SUB_FFF, "sub_fff", Tf, Tf, Tf, T0)
+OPCODE(71, OP_MUL_FFF, "mul_fff", Tf, Tf, Tf, T0)
+OPCODE(72, OP_DIV_FFF, "div_fff", Tf, Tf, Tf, T0)
+OPCODE(73, OP_POW_FFF, "pow_fff", Tf, Tf, Tf, T0)
+OPCODE(74, OP_MOD_FFF, "mod_fff", Tf, Tf, Tf, T0)
+OPCODE(75, OP_FLOORDIV_FFF, "floordiv_fff", Tf, Tf, Tf, T0)
+OPCODE(76, OP_SQRT_FF, "sqrt_ff", Tf, Tf, T0, T0)
+OPCODE(77, OP_WHERE_FBFF, "where_fbff", Tf, Tb, Tf, Tf)
+
+OPCODE(78, OP_FUNC_FFN, "func_ffn", Tf, Tf, Tn, T0)
+OPCODE(79, OP_FUNC_FFFN, "func_fffn", Tf, Tf, Tf, Tn)
+
+OPCODE(80, OP_CAST_DI, "cast_di", Td, Ti, T0, T0)
+OPCODE(81, OP_CAST_DL, "cast_dl", Td, Tl, T0, T0)
+OPCODE(82, OP_CAST_DF, "cast_df", Td, Tf, T0, T0)
+OPCODE(83, OP_COPY_DD, "copy_dd", Td, Td, T0, T0)
+OPCODE(84, OP_ONES_LIKE_DD, "ones_like_dd", Td, T0, T0, T0)
+OPCODE(85, OP_NEG_DD, "neg_dd", Td, Td, T0, T0)
+OPCODE(86, OP_ADD_DDD, "add_ddd", Td, Td, Td, T0)
+OPCODE(87, OP_SUB_DDD, "sub_ddd", Td, Td, Td, T0)
+OPCODE(88, OP_MUL_DDD, "mul_ddd", Td, Td, Td, T0)
+OPCODE(89, OP_DIV_DDD, "div_ddd", Td, Td, Td, T0)
+OPCODE(90, OP_POW_DDD, "pow_ddd", Td, Td, Td, T0)
+OPCODE(91, OP_MOD_DDD, "mod_ddd", Td, Td, Td, T0)
+OPCODE(92, OP_FLOORDIV_DDD, "floordiv_ddd", Td, Td, Td, T0)
+
+OPCODE(93, OP_SQRT_DD, "sqrt_dd", Td, Td, T0, T0)
+OPCODE(94, OP_WHERE_DBDD, "where_dbdd", Td, Tb, Td, Td)
+OPCODE(95, OP_FUNC_DDN, "func_ddn", Td, Td, Tn, T0)
+OPCODE(96, OP_FUNC_DDDN, "func_dddn", Td, Td, Td, Tn)
+
+OPCODE(97, OP_EQ_BCC, "eq_bcc", Tb, Tc, Tc, T0)
+OPCODE(98, OP_NE_BCC, "ne_bcc", Tb, Tc, Tc, T0)
+
+OPCODE(99, OP_CAST_CI, "cast_ci", Tc, Ti, T0, T0)
+OPCODE(100, OP_CAST_CL, "cast_cl", Tc, Tl, T0, T0)
+OPCODE(101, OP_CAST_CF, "cast_cf", Tc, Tf, T0, T0)
+OPCODE(102, OP_CAST_CD, "cast_cd", Tc, Td, T0, T0)
+OPCODE(103, OP_ONES_LIKE_CC, "ones_like_cc", Tc, T0, T0, T0)
+OPCODE(104, OP_COPY_CC, "copy_cc", Tc, Tc, T0, T0)
+OPCODE(105, OP_NEG_CC, "neg_cc", Tc, Tc, T0, T0)
+OPCODE(106, OP_ADD_CCC, "add_ccc", Tc, Tc, Tc, T0)
+OPCODE(107, OP_SUB_CCC, "sub_ccc", Tc, Tc, Tc, T0)
+OPCODE(108, OP_MUL_CCC, "mul_ccc", Tc, Tc, Tc, T0)
+OPCODE(109, OP_DIV_CCC, "div_ccc", Tc, Tc, Tc, T0)
+OPCODE(110, OP_WHERE_CBCC, "where_cbcc", Tc, Tb, Tc, Tc)
+OPCODE(111, OP_FUNC_CCN, "func_ccn", Tc, Tc, Tn, T0)
+OPCODE(112, OP_FUNC_CCCN, "func_cccn", Tc, Tc, Tc, Tn)
+
+OPCODE(113, OP_REAL_DC, "real_dc", Td, Tc, T0, T0)
+OPCODE(114, OP_IMAG_DC, "imag_dc", Td, Tc, T0, T0)
+OPCODE(115, OP_COMPLEX_CDD, "complex_cdd", Tc, Td, Td, T0)
+
+OPCODE(116, OP_COPY_SS, "copy_ss", Ts, Ts, T0, T0)
+
+OPCODE(117, OP_WHERE_BBBB, "where_bbbb", Tb, Tb, Tb, Tb)
+
+OPCODE(118, OP_CONTAINS_BSS, "contains_bss", Tb, Ts, Ts, T0)
+//Boolean outputs
+OPCODE(119, OP_FUNC_BDN, "func_bdn", Tb, Td, Tn, T0)
+OPCODE(120, OP_FUNC_BFN, "func_bfn", Tb, Tf, Tn, T0)
+OPCODE(121, OP_FUNC_BCN, "func_bcn", Tb, Tc, Tn, T0)
+//Integer funcs
+OPCODE(122, OP_FUNC_IIN, "func_iin", Ti, Ti, Tn, T0)
+OPCODE(123, OP_FUNC_LLN, "func_lln", Tl, Tl, Tn, T0)
+
+// Reductions always have to be at the end - parts of the code
+// use > OP_REDUCTION to decide whether operation is a reduction
+OPCODE(124, OP_REDUCTION, NULL, T0, T0, T0, T0)
+
+/* Last argument in a reduction is the axis of the array the
+   reduction should be applied along. */
+
+OPCODE(125, OP_SUM_IIN, "sum_iin", Ti, Ti, Tn, T0)
+OPCODE(126, OP_SUM_LLN, "sum_lln", Tl, Tl, Tn, T0)
+OPCODE(127, OP_SUM_FFN, "sum_ffn", Tf, Tf, Tn, T0)
+OPCODE(128, OP_SUM_DDN, "sum_ddn", Td, Td, Tn, T0)
+OPCODE(129, OP_SUM_CCN, "sum_ccn", Tc, Tc, Tn, T0)
+
+OPCODE(130, OP_PROD, NULL, T0, T0, T0, T0)
+OPCODE(131, OP_PROD_IIN, "prod_iin", Ti, Ti, Tn, T0)
+OPCODE(132, OP_PROD_LLN, "prod_lln", Tl, Tl, Tn, T0)
+OPCODE(133, OP_PROD_FFN, "prod_ffn", Tf, Tf, Tn, T0)
+OPCODE(134, OP_PROD_DDN, "prod_ddn", Td, Td, Tn, T0)
+OPCODE(135, OP_PROD_CCN, "prod_ccn", Tc, Tc, Tn, T0)
+
+OPCODE(136, OP_MIN, NULL, T0, T0, T0, T0)
+OPCODE(137, OP_MIN_IIN, "min_iin", Ti, Ti, Tn, T0)
+OPCODE(138, OP_MIN_LLN, "min_lln", Tl, Tl, Tn, T0)
+OPCODE(139, OP_MIN_FFN, "min_ffn", Tf, Tf, Tn, T0)
+OPCODE(140, OP_MIN_DDN, "min_ddn", Td, Td, Tn, T0)
+
+OPCODE(141, OP_MAX, NULL, T0, T0, T0, T0)
+OPCODE(142, OP_MAX_IIN, "max_iin", Ti, Ti, Tn, T0)
+OPCODE(143, OP_MAX_LLN, "max_lln", Tl, Tl, Tn, T0)
+OPCODE(144, OP_MAX_FFN, "max_ffn", Tf, Tf, Tn, T0)
+OPCODE(145, OP_MAX_DDN, "max_ddn", Td, Td, Tn, T0)
+
+/*
+When we get to 255, will maybe have to change code again
+(change latin_1 encoding in necompiler.py, use something
+other than unsigned char for OPCODE table)
+*/
+/* Should be the last opcode */
+OPCODE(146, OP_END, NULL, T0, T0, T0, T0)

numexpr/str-two-way.hpp

@@ -0,0 +1,435 @@
+/* Byte-wise substring search, using the Two-Way algorithm.
+ * Copyright (C) 2008, 2010 Eric Blake
+ * Permission to use, copy, modify, and distribute this software
+ * is freely granted, provided that this notice is preserved.
+ */
+
+
+/* Before including this file, you need to include <string.h>, and define:
+     RETURN_TYPE		A macro that expands to the return type.
+     AVAILABLE(h, h_l, j, n_l)	A macro that returns nonzero if there are
+				at least N_L bytes left starting at
+				H[J].  H is 'unsigned char *', H_L, J,
+				and N_L are 'size_t'; H_L is an
+				lvalue.  For NUL-terminated searches,
+				H_L can be modified each iteration to
+				avoid having to compute the end of H
+				up front.
+
+  For case-insensitivity, you may optionally define:
+     CMP_FUNC(p1, p2, l)	A macro that returns 0 iff the first L
+				characters of P1 and P2 are equal.
+     CANON_ELEMENT(c)		A macro that canonicalizes an element
+				right after it has been fetched from
+				one of the two strings.  The argument
+				is an 'unsigned char'; the result must
+				be an 'unsigned char' as well.
+
+  This file undefines the macros documented above, and defines
+  LONG_NEEDLE_THRESHOLD.
+*/
+
+#include <limits.h>
+
+/*
+  Python 2.7 (the only Python 2.x version supported as of now and until 2020)
+  is built on windows with Visual Studio 2008 C compiler. That dictates that
+  the compiler which must be used by authors of third party Python modules.
+  See https://mail.python.org/pipermail/distutils-sig/2014-September/024885.html
+
+  Unfortunately this version of Visual Studio doesn't claim to be C99 compatible
+  and in particular it lacks the stdint.h header. So we have to replace it with
+  a public domain version.
+
+  Visual Studio 2010 and later have stdint.h.
+*/
+
+#ifdef _MSC_VER
+  #if _MSC_VER <= 1500
+    #include "win32/stdint.h"
+  #endif
+#else
+  #include <stdint.h>
+#endif
+
+/* We use the Two-Way string matching algorithm, which guarantees
+   linear complexity with constant space.  Additionally, for long
+   needles, we also use a bad character shift table similar to the
+   Boyer-Moore algorithm to achieve improved (potentially sub-linear)
+   performance.
+
+   See http://www-igm.univ-mlv.fr/~lecroq/string/node26.html#SECTION00260
+   and http://en.wikipedia.org/wiki/Boyer-Moore_string_search_algorithm
+*/
+
+/* Point at which computing a bad-byte shift table is likely to be
+   worthwhile.  Small needles should not compute a table, since it
+   adds (1 << CHAR_BIT) + NEEDLE_LEN computations of preparation for a
+   speedup no greater than a factor of NEEDLE_LEN.  The larger the
+   needle, the better the potential performance gain.  On the other
+   hand, on non-POSIX systems with CHAR_BIT larger than eight, the
+   memory required for the table is prohibitive.  */
+#if CHAR_BIT < 10
+# define LONG_NEEDLE_THRESHOLD 32U
+#else
+# define LONG_NEEDLE_THRESHOLD SIZE_MAX
+#endif
+
+#define MAX(a, b) ((a < b) ? (b) : (a))
+
+#ifndef CANON_ELEMENT
+# define CANON_ELEMENT(c) c
+#endif
+#ifndef CMP_FUNC
+# define CMP_FUNC memcmp
+#endif
+
+/* Perform a critical factorization of NEEDLE, of length NEEDLE_LEN.
+   Return the index of the first byte in the right half, and set
+   *PERIOD to the global period of the right half.
+
+   The global period of a string is the smallest index (possibly its
+   length) at which all remaining bytes in the string are repetitions
+   of the prefix (the last repetition may be a subset of the prefix).
+
+   When NEEDLE is factored into two halves, a local period is the
+   length of the smallest word that shares a suffix with the left half
+   and shares a prefix with the right half.  All factorizations of a
+   non-empty NEEDLE have a local period of at least 1 and no greater
+   than NEEDLE_LEN.
+
+   A critical factorization has the property that the local period
+   equals the global period.  All strings have at least one critical
+   factorization with the left half smaller than the global period.
+
+   Given an ordered alphabet, a critical factorization can be computed
+   in linear time, with 2 * NEEDLE_LEN comparisons, by computing the
+   larger of two ordered maximal suffixes.  The ordered maximal
+   suffixes are determined by lexicographic comparison of
+   periodicity.  */
+static size_t
+critical_factorization (const unsigned char *needle, size_t needle_len,
+			size_t *period)
+{
+  /* Index of last byte of left half, or SIZE_MAX.  */
+  size_t max_suffix, max_suffix_rev;
+  size_t j; /* Index into NEEDLE for current candidate suffix.  */
+  size_t k; /* Offset into current period.  */
+  size_t p; /* Intermediate period.  */
+  unsigned char a, b; /* Current comparison bytes.  */
+
+  /* Invariants:
+     0 <= j < NEEDLE_LEN - 1
+     -1 <= max_suffix{,_rev} < j (treating SIZE_MAX as if it were signed)
+     min(max_suffix, max_suffix_rev) < global period of NEEDLE
+     1 <= p <= global period of NEEDLE
+     p == global period of the substring NEEDLE[max_suffix{,_rev}+1...j]
+     1 <= k <= p
+  */
+
+  /* Perform lexicographic search.  */
+  max_suffix = SIZE_MAX;
+  j = 0;
+  k = p = 1;
+  while (j + k < needle_len)
+    {
+      a = CANON_ELEMENT (needle[j + k]);
+      b = CANON_ELEMENT (needle[(size_t)(max_suffix + k)]);
+      if (a < b)
+	{
+	  /* Suffix is smaller, period is entire prefix so far.  */
+	  j += k;
+	  k = 1;
+	  p = j - max_suffix;
+	}
+      else if (a == b)
+	{
+	  /* Advance through repetition of the current period.  */
+	  if (k != p)
+	    ++k;
+	  else
+	    {
+	      j += p;
+	      k = 1;
+	    }
+	}
+      else /* b < a */
+	{
+	  /* Suffix is larger, start over from current location.  */
+	  max_suffix = j++;
+	  k = p = 1;
+	}
+    }
+  *period = p;
+
+  /* Perform reverse lexicographic search.  */
+  max_suffix_rev = SIZE_MAX;
+  j = 0;
+  k = p = 1;
+  while (j + k < needle_len)
+    {
+      a = CANON_ELEMENT (needle[j + k]);
+      b = CANON_ELEMENT (needle[max_suffix_rev + k]);
+      if (b < a)
+	{
+	  /* Suffix is smaller, period is entire prefix so far.  */
+	  j += k;
+	  k = 1;
+	  p = j - max_suffix_rev;
+	}
+      else if (a == b)
+	{
+	  /* Advance through repetition of the current period.  */
+	  if (k != p)
+	    ++k;
+	  else
+	    {
+	      j += p;
+	      k = 1;
+	    }
+	}
+      else /* a < b */
+	{
+	  /* Suffix is larger, start over from current location.  */
+	  max_suffix_rev = j++;
+	  k = p = 1;
+	}
+    }
+
+  /* Choose the longer suffix.  Return the first byte of the right
+     half, rather than the last byte of the left half.  */
+  if (max_suffix_rev + 1 < max_suffix + 1)
+    return max_suffix + 1;
+  *period = p;
+  return max_suffix_rev + 1;
+}
+
+/* Return the first location of non-empty NEEDLE within HAYSTACK, or
+   NULL.  HAYSTACK_LEN is the minimum known length of HAYSTACK.  This
+   method is optimized for NEEDLE_LEN < LONG_NEEDLE_THRESHOLD.
+   Performance is guaranteed to be linear, with an initialization cost
+   of 2 * NEEDLE_LEN comparisons.
+
+   If AVAILABLE does not modify HAYSTACK_LEN (as in memmem), then at
+   most 2 * HAYSTACK_LEN - NEEDLE_LEN comparisons occur in searching.
+   If AVAILABLE modifies HAYSTACK_LEN (as in strstr), then at most 3 *
+   HAYSTACK_LEN - NEEDLE_LEN comparisons occur in searching.  */
+static RETURN_TYPE
+two_way_short_needle (const unsigned char *haystack, size_t haystack_len,
+		      const unsigned char *needle, size_t needle_len)
+{
+  size_t i; /* Index into current byte of NEEDLE.  */
+  size_t j; /* Index into current window of HAYSTACK.  */
+  size_t period; /* The period of the right half of needle.  */
+  size_t suffix; /* The index of the right half of needle.  */
+
+  /* Factor the needle into two halves, such that the left half is
+     smaller than the global period, and the right half is
+     periodic (with a period as large as NEEDLE_LEN - suffix).  */
+  suffix = critical_factorization (needle, needle_len, &period);
+
+  /* Perform the search.  Each iteration compares the right half
+     first.  */
+  if (CMP_FUNC (needle, needle + period, suffix) == 0)
+    {
+      /* Entire needle is periodic; a mismatch can only advance by the
+	 period, so use memory to avoid rescanning known occurrences
+	 of the period.  */
+      size_t memory = 0;
+      j = 0;
+      while (AVAILABLE (haystack, haystack_len, j, needle_len))
+	{
+	  /* Scan for matches in right half.  */
+	  i = MAX (suffix, memory);
+	  while (i < needle_len && (CANON_ELEMENT (needle[i])
+				    == CANON_ELEMENT (haystack[i + j])))
+	    ++i;
+	  if (needle_len <= i)
+	    {
+	      /* Scan for matches in left half.  */
+	      i = suffix - 1;
+	      while (memory < i + 1 && (CANON_ELEMENT (needle[i])
+					== CANON_ELEMENT (haystack[i + j])))
+		--i;
+	      if (i + 1 < memory + 1)
+		return (RETURN_TYPE) (haystack + j);
+	      /* No match, so remember how many repetitions of period
+		 on the right half were scanned.  */
+	      j += period;
+	      memory = needle_len - period;
+	    }
+	  else
+	    {
+	      j += i - suffix + 1;
+	      memory = 0;
+	    }
+	}
+    }
+  else
+    {
+      /* The two halves of needle are distinct; no extra memory is
+	 required, and any mismatch results in a maximal shift.  */
+      period = MAX (suffix, needle_len - suffix) + 1;
+      j = 0;
+      while (AVAILABLE (haystack, haystack_len, j, needle_len))
+	{
+	  /* Scan for matches in right half.  */
+	  i = suffix;
+	  while (i < needle_len && (CANON_ELEMENT (needle[i])
+				    == CANON_ELEMENT (haystack[i + j])))
+	    ++i;
+	  if (needle_len <= i)
+	    {
+	      /* Scan for matches in left half.  */
+	      i = suffix - 1;
+	      while (i != SIZE_MAX && (CANON_ELEMENT (needle[i])
+				       == CANON_ELEMENT (haystack[i + j])))
+		--i;
+	      if (i == SIZE_MAX)
+		return (RETURN_TYPE) (haystack + j);
+	      j += period;
+	    }
+	  else
+	    j += i - suffix + 1;
+	}
+    }
+  return NULL;
+}
+
+/* Return the first location of non-empty NEEDLE within HAYSTACK, or
+   NULL.  HAYSTACK_LEN is the minimum known length of HAYSTACK.  This
+   method is optimized for LONG_NEEDLE_THRESHOLD <= NEEDLE_LEN.
+   Performance is guaranteed to be linear, with an initialization cost
+   of 3 * NEEDLE_LEN + (1 << CHAR_BIT) operations.
+
+   If AVAILABLE does not modify HAYSTACK_LEN (as in memmem), then at
+   most 2 * HAYSTACK_LEN - NEEDLE_LEN comparisons occur in searching,
+   and sublinear performance O(HAYSTACK_LEN / NEEDLE_LEN) is possible.
+   If AVAILABLE modifies HAYSTACK_LEN (as in strstr), then at most 3 *
+   HAYSTACK_LEN - NEEDLE_LEN comparisons occur in searching, and
+   sublinear performance is not possible.  */
+static RETURN_TYPE
+two_way_long_needle (const unsigned char *haystack, size_t haystack_len,
+		     const unsigned char *needle, size_t needle_len)
+{
+  size_t i; /* Index into current byte of NEEDLE.  */
+  size_t j; /* Index into current window of HAYSTACK.  */
+  size_t period; /* The period of the right half of needle.  */
+  size_t suffix; /* The index of the right half of needle.  */
+  size_t shift_table[1U << CHAR_BIT]; /* See below.  */
+
+  /* Factor the needle into two halves, such that the left half is
+     smaller than the global period, and the right half is
+     periodic (with a period as large as NEEDLE_LEN - suffix).  */
+  suffix = critical_factorization (needle, needle_len, &period);
+
+  /* Populate shift_table.  For each possible byte value c,
+     shift_table[c] is the distance from the last occurrence of c to
+     the end of NEEDLE, or NEEDLE_LEN if c is absent from the NEEDLE.
+     shift_table[NEEDLE[NEEDLE_LEN - 1]] contains the only 0.  */
+  for (i = 0; i < 1U << CHAR_BIT; i++)
+    shift_table[i] = needle_len;
+  for (i = 0; i < needle_len; i++)
+    shift_table[CANON_ELEMENT (needle[i])] = needle_len - i - 1;
+
+  /* Perform the search.  Each iteration compares the right half
+     first.  */
+  if (CMP_FUNC (needle, needle + period, suffix) == 0)
+    {
+      /* Entire needle is periodic; a mismatch can only advance by the
+	 period, so use memory to avoid rescanning known occurrences
+	 of the period.  */
+      size_t memory = 0;
+      size_t shift;
+      j = 0;
+      while (AVAILABLE (haystack, haystack_len, j, needle_len))
+	{
+	  /* Check the last byte first; if it does not match, then
+	     shift to the next possible match location.  */
+	  shift = shift_table[CANON_ELEMENT (haystack[j + needle_len - 1])];
+	  if (0 < shift)
+	    {
+	      if (memory && shift < period)
+		{
+		  /* Since needle is periodic, but the last period has
+		     a byte out of place, there can be no match until
+		     after the mismatch.  */
+		  shift = needle_len - period;
+		}
+	      memory = 0;
+	      j += shift;
+	      continue;
+	    }
+	  /* Scan for matches in right half.  The last byte has
+	     already been matched, by virtue of the shift table.  */
+	  i = MAX (suffix, memory);
+	  while (i < needle_len - 1 && (CANON_ELEMENT (needle[i])
+					== CANON_ELEMENT (haystack[i + j])))
+	    ++i;
+	  if (needle_len - 1 <= i)
+	    {
+	      /* Scan for matches in left half.  */
+	      i = suffix - 1;
+	      while (memory < i + 1 && (CANON_ELEMENT (needle[i])
+					== CANON_ELEMENT (haystack[i + j])))
+		--i;
+	      if (i + 1 < memory + 1)
+		return (RETURN_TYPE) (haystack + j);
+	      /* No match, so remember how many repetitions of period
+		 on the right half were scanned.  */
+	      j += period;
+	      memory = needle_len - period;
+	    }
+	  else
+	    {
+	      j += i - suffix + 1;
+	      memory = 0;
+	    }
+	}
+    }
+  else
+    {
+      /* The two halves of needle are distinct; no extra memory is
+	 required, and any mismatch results in a maximal shift.  */
+      size_t shift;
+      period = MAX (suffix, needle_len - suffix) + 1;
+      j = 0;
+      while (AVAILABLE (haystack, haystack_len, j, needle_len))
+	{
+	  /* Check the last byte first; if it does not match, then
+	     shift to the next possible match location.  */
+	  shift = shift_table[CANON_ELEMENT (haystack[j + needle_len - 1])];
+	  if (0 < shift)
+	    {
+	      j += shift;
+	      continue;
+	    }
+	  /* Scan for matches in right half.  The last byte has
+	     already been matched, by virtue of the shift table.  */
+	  i = suffix;
+	  while (i < needle_len - 1 && (CANON_ELEMENT (needle[i])
+					== CANON_ELEMENT (haystack[i + j])))
+	    ++i;
+	  if (needle_len - 1 <= i)
+	    {
+	      /* Scan for matches in left half.  */
+	      i = suffix - 1;
+	      while (i != SIZE_MAX && (CANON_ELEMENT (needle[i])
+				       == CANON_ELEMENT (haystack[i + j])))
+		--i;
+	      if (i == SIZE_MAX)
+		return (RETURN_TYPE) (haystack + j);
+	      j += period;
+	    }
+	  else
+	    j += i - suffix + 1;
+	}
+    }
+  return NULL;
+}
+
+#undef AVAILABLE
+#undef CANON_ELEMENT
+#undef CMP_FUNC
+#undef MAX
+#undef RETURN_TYPE

numexpr/tests/__init__.py

@@ -0,0 +1,14 @@
+###################################################################
+#  Numexpr - Fast numerical array expression evaluator for NumPy.
+#
+#      License: MIT
+#      Author:  See AUTHORS.txt
+#
+#  See LICENSE.txt and LICENSES/*.txt for details about copyright and
+#  rights to use.
+####################################################################
+
+from numexpr.tests.test_numexpr import print_versions, test
+
+if __name__ == '__main__':
+    test()

numexpr/tests/conftest.py

@@ -0,0 +1,21 @@
+###################################################################
+#  Numexpr - Fast numerical array expression evaluator for NumPy.
+#
+#      License: MIT
+#      Author:  See AUTHORS.txt
+#
+#  See LICENSE.txt and LICENSES/*.txt for details about copyright and
+#  rights to use.
+####################################################################
+
+import pytest
+
+import numexpr
+
+
+def pytest_configure(config):
+    config.addinivalue_line(
+        "markers", "thread_unsafe: mark test as unsafe for parallel execution"
+    )
+    print("")
+    numexpr.print_versions()

numexpr/utils.py

@@ -0,0 +1,311 @@
+###################################################################
+#  Numexpr - Fast numerical array expression evaluator for NumPy.
+#
+#      License: MIT
+#      Author:  See AUTHORS.txt
+#
+#  See LICENSE.txt and LICENSES/*.txt for details about copyright and
+#  rights to use.
+####################################################################
+
+import logging
+
+log = logging.getLogger(__name__)
+
+import contextvars
+import os
+import subprocess
+
+from numexpr import use_vml
+from numexpr.interpreter import MAX_THREADS, _get_num_threads, _set_num_threads
+
+from . import version
+
+if use_vml:
+    from numexpr.interpreter import (_get_vml_num_threads, _get_vml_version,
+                                     _set_vml_accuracy_mode,
+                                     _set_vml_num_threads)
+
+
+def get_vml_version():
+    """
+    Get the VML/MKL library version.
+    """
+    if use_vml:
+        return _get_vml_version()
+    else:
+        return None
+
+
+def set_vml_accuracy_mode(mode):
+    """
+    Set the accuracy mode for VML operations.
+
+    The `mode` parameter can take the values:
+    - 'high': high accuracy mode (HA), <1 least significant bit
+    - 'low': low accuracy mode (LA), typically 1-2 least significant bits
+    - 'fast': enhanced performance mode (EP)
+    - None: mode settings are ignored
+
+    This call is equivalent to the `vmlSetMode()` in the VML library.
+    See:
+
+    http://www.intel.com/software/products/mkl/docs/webhelp/vml/vml_DataTypesAccuracyModes.html
+
+    for more info on the accuracy modes.
+
+    Returns old accuracy settings.
+    """
+    if use_vml:
+        acc_dict = {None: 0, 'low': 1, 'high': 2, 'fast': 3}
+        acc_reverse_dict = {1: 'low', 2: 'high', 3: 'fast'}
+        if mode not in list(acc_dict.keys()):
+            raise ValueError(
+                "mode argument must be one of: None, 'high', 'low', 'fast'")
+        retval = _set_vml_accuracy_mode(acc_dict.get(mode, 0))
+        return acc_reverse_dict.get(retval)
+    else:
+        return None
+
+
+def set_vml_num_threads(nthreads):
+    """
+    Suggests a maximum number of threads to be used in VML operations.
+
+    This function is equivalent to the call
+    `mkl_domain_set_num_threads(nthreads, MKL_DOMAIN_VML)` in the MKL
+    library.  See:
+
+    http://www.intel.com/software/products/mkl/docs/webhelp/support/functn_mkl_domain_set_num_threads.html
+
+    for more info about it.
+    """
+    if use_vml:
+        _set_vml_num_threads(nthreads)
+    pass
+
+def get_vml_num_threads():
+    """
+    Gets the maximum number of threads to be used in VML operations.
+
+    This function is equivalent to the call
+    `mkl_domain_get_max_threads (MKL_DOMAIN_VML)` in the MKL
+    library.  See:
+
+    http://software.intel.com/en-us/node/522118
+
+    for more info about it.
+    """
+    if use_vml:
+        return _get_vml_num_threads()
+    return None
+
+def set_num_threads(nthreads):
+    """
+    Sets a number of threads to be used in operations.
+
+    DEPRECATED: returns the previous setting for the number of threads.
+
+    During initialization time NumExpr sets this number to the number
+    of detected cores in the system (see `detect_number_of_cores()`).
+    """
+    old_nthreads = _set_num_threads(nthreads)
+    return old_nthreads
+
+def get_num_threads():
+    """
+    Gets the number of threads currently in use for operations.
+    """
+    return _get_num_threads()
+
+def _init_num_threads():
+    """
+    Detects the environment variable 'NUMEXPR_MAX_THREADS' to set the threadpool
+    size, and if necessary the slightly redundant 'NUMEXPR_NUM_THREADS' or
+    'OMP_NUM_THREADS' env vars to set the initial number of threads used by
+    the virtual machine.
+    """
+    # Any platform-specific short-circuits
+    if 'sparc' in version.platform_machine:
+        log.warning('The number of threads have been set to 1 because problems related '
+                  'to threading have been reported on some sparc machine. '
+                  'The number of threads can be changed using the "set_num_threads" '
+                  'function.')
+        set_num_threads(1)
+        return 1
+
+    env_configured = False
+    n_cores = detect_number_of_cores()
+    if ('NUMEXPR_MAX_THREADS' in os.environ and os.environ['NUMEXPR_MAX_THREADS'] != '' or
+        'OMP_NUM_THREADS' in os.environ and os.environ['OMP_NUM_THREADS'] != ''):
+        # The user has configured NumExpr in the expected way, so suppress logs.
+        env_configured = True
+        n_cores = MAX_THREADS
+    else:
+        # The use has not set 'NUMEXPR_MAX_THREADS', so likely they have not
+        # configured NumExpr as desired, so we emit info logs.
+        if n_cores > MAX_THREADS:
+            log.info('Note: detected %d virtual cores but NumExpr set to maximum of %d, check "NUMEXPR_MAX_THREADS" environment variable.'%(n_cores, MAX_THREADS))
+        if n_cores > 16:
+            # Back in 2019, 8 threads would be considered safe for performance. We are in 2024 now, so adjusting.
+            log.info('Note: NumExpr detected %d cores but "NUMEXPR_MAX_THREADS" not set, so enforcing safe limit of 16.'%n_cores)
+            n_cores = 16
+
+    # Now we check for 'NUMEXPR_NUM_THREADS' or 'OMP_NUM_THREADS' to set the
+    # actual number of threads used.
+    if 'NUMEXPR_NUM_THREADS' in os.environ and os.environ['NUMEXPR_NUM_THREADS'] != '':
+        requested_threads = int(os.environ['NUMEXPR_NUM_THREADS'])
+    elif 'OMP_NUM_THREADS' in os.environ and os.environ['OMP_NUM_THREADS'] != '':
+        # Empty string is commonly used to unset the variable
+        requested_threads = int(os.environ['OMP_NUM_THREADS'])
+    else:
+        requested_threads = n_cores
+        if not env_configured:
+            log.info('NumExpr defaulting to %d threads.'%n_cores)
+
+    # The C-extension function performs its own checks against `MAX_THREADS`
+    set_num_threads(requested_threads)
+    return requested_threads
+
+
+def detect_number_of_cores():
+    """
+    Detects the number of cores on a system. Cribbed from pp.
+    """
+    # Linux, Unix and MacOS:
+    if hasattr(os, "sysconf"):
+        if "SC_NPROCESSORS_ONLN" in os.sysconf_names:
+            # Linux & Unix:
+            ncpus = os.sysconf("SC_NPROCESSORS_ONLN")
+            if isinstance(ncpus, int) and ncpus > 0:
+                return ncpus
+        else:  # OSX:
+            return int(subprocess.check_output(["sysctl", "-n", "hw.ncpu"]))
+    # Windows:
+    try:
+        ncpus = int(os.environ.get("NUMBER_OF_PROCESSORS", ""))
+        if ncpus > 0:
+            return ncpus
+    except ValueError:
+        pass
+    return 1  # Default
+
+
+def detect_number_of_threads():
+    """
+    DEPRECATED: use `_init_num_threads` instead.
+    If this is modified, please update the note in: https://github.com/pydata/numexpr/wiki/Numexpr-Users-Guide
+    """
+    log.warning('Deprecated, use `init_num_threads` instead.')
+    try:
+        nthreads = int(os.environ.get('NUMEXPR_NUM_THREADS', ''))
+    except ValueError:
+        try:
+            nthreads = int(os.environ.get('OMP_NUM_THREADS', ''))
+        except ValueError:
+            nthreads = detect_number_of_cores()
+
+    # Check that we don't surpass the MAX_THREADS in interpreter.cpp
+    if nthreads > MAX_THREADS:
+        nthreads = MAX_THREADS
+    return nthreads
+
+
+class CacheDict(dict):
+    """
+    A dictionary that prevents itself from growing too much.
+    """
+
+    def __init__(self, maxentries):
+        self.maxentries = maxentries
+        super(CacheDict, self).__init__(self)
+
+    def __setitem__(self, key, value):
+        # Protection against growing the cache too much
+        if len(self) > self.maxentries:
+            # Remove a 10% of (arbitrary) elements from the cache
+            entries_to_remove = self.maxentries // 10
+            for k in list(self.keys())[:entries_to_remove]:
+                super(CacheDict, self).__delitem__(k)
+        super(CacheDict, self).__setitem__(key, value)
+
+
+class ContextDict:
+    """
+    A context aware version dictionary
+    """
+    def __init__(self):
+        self._context_data = contextvars.ContextVar('context_data', default={})
+
+    def set(self, key=None, value=None, **kwargs):
+        data = self._context_data.get().copy()
+
+        if key is not None:
+            data[key] = value
+
+        for k, v in kwargs.items():
+            data[k] = v
+
+        self._context_data.set(data)
+
+    def get(self, key, default=None):
+        data = self._context_data.get()
+        return data.get(key, default)
+
+    def delete(self, key):
+        data = self._context_data.get().copy()
+        if key in data:
+            del data[key]
+        self._context_data.set(data)
+
+    def clear(self):
+        self._context_data.set({})
+
+    def all(self):
+        return self._context_data.get()
+
+    def update(self, *args, **kwargs):
+        data = self._context_data.get().copy()
+
+        if args:
+            if len(args) > 1:
+                raise TypeError(f"update() takes at most 1 positional argument ({len(args)} given)")
+            other = args[0]
+            if isinstance(other, dict):
+                data.update(other)
+            else:
+                for k, v in other:
+                    data[k] = v
+
+        data.update(kwargs)
+        self._context_data.set(data)
+
+    def keys(self):
+        return self._context_data.get().keys()
+
+    def values(self):
+        return self._context_data.get().values()
+
+    def items(self):
+        return self._context_data.get().items()
+
+    def __getitem__(self, key):
+        return self.get(key)
+
+    def __setitem__(self, key, value):
+        self.set(key, value)
+
+    def __delitem__(self, key):
+        self.delete(key)
+
+    def __contains__(self, key):
+        return key in self._context_data.get()
+
+    def __len__(self):
+        return len(self._context_data.get())
+
+    def __iter__(self):
+        return iter(self._context_data.get())
+
+    def __repr__(self):
+        return repr(self._context_data.get())

numexpr/win32/pthread.c

@@ -0,0 +1,218 @@
+/*
+ * Code for simulating pthreads API on Windows.  This is Git-specific,
+ * but it is enough for Numexpr needs too.
+ *
+ * Copyright (C) 2009 Andrzej K. Haczewski <ahaczewski@gmail.com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ *
+ * DISCLAIMER: The implementation is Git-specific, it is subset of original
+ * Pthreads API, without lots of other features that Git doesn't use.
+ * Git also makes sure that the passed arguments are valid, so there's
+ * no need for double-checking.
+ */
+
+#include "pthread.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <errno.h>
+#include <limits.h>
+#include <process.h>
+
+
+void die(const char *err, ...)
+{
+    printf("%s", err);
+    exit(-1);
+}
+
+static unsigned __stdcall win32_start_routine(void *arg)
+{
+    pthread_t *thread = arg;
+    thread->arg = thread->start_routine(thread->arg);
+    return 0;
+}
+
+int pthread_create(pthread_t *thread, const void *unused,
+           void *(*start_routine)(void*), void *arg)
+{
+    thread->arg = arg;
+    thread->start_routine = start_routine;
+    thread->handle = (HANDLE)
+        _beginthreadex(NULL, 0, win32_start_routine, thread, 0, NULL);
+
+    if (!thread->handle)
+        return errno;
+    else
+        return 0;
+}
+
+int win32_pthread_join(pthread_t *thread, void **value_ptr)
+{
+    DWORD result = WaitForSingleObject(thread->handle, INFINITE);
+    switch (result) {
+        case WAIT_OBJECT_0:
+            if (value_ptr)
+                *value_ptr = thread->arg;
+            return 0;
+        case WAIT_ABANDONED:
+            return EINVAL;
+        default:
+            return GetLastError();
+    }
+}
+
+int pthread_cond_init(pthread_cond_t *cond, const void *unused)
+{
+    cond->waiters = 0;
+    cond->was_broadcast = 0;
+    InitializeCriticalSection(&cond->waiters_lock);
+
+    cond->sema = CreateSemaphore(NULL, 0, LONG_MAX, NULL);
+    if (!cond->sema)
+        die("CreateSemaphore() failed");
+
+    cond->continue_broadcast = CreateEvent(NULL,    /* security */
+                FALSE,            /* auto-reset */
+                FALSE,            /* not signaled */
+                NULL);            /* name */
+    if (!cond->continue_broadcast)
+        die("CreateEvent() failed");
+
+    return 0;
+}
+
+int pthread_cond_destroy(pthread_cond_t *cond)
+{
+    CloseHandle(cond->sema);
+    CloseHandle(cond->continue_broadcast);
+    DeleteCriticalSection(&cond->waiters_lock);
+    return 0;
+}
+
+int pthread_cond_wait(pthread_cond_t *cond, CRITICAL_SECTION *mutex)
+{
+    int last_waiter;
+
+    EnterCriticalSection(&cond->waiters_lock);
+    cond->waiters++;
+    LeaveCriticalSection(&cond->waiters_lock);
+
+    /*
+     * Unlock external mutex and wait for signal.
+     * NOTE: we've held mutex locked long enough to increment
+     * waiters count above, so there's no problem with
+     * leaving mutex unlocked before we wait on semaphore.
+     */
+    LeaveCriticalSection(mutex);
+
+    /* let's wait - ignore return value */
+    WaitForSingleObject(cond->sema, INFINITE);
+
+    /*
+     * Decrease waiters count. If we are the last waiter, then we must
+     * notify the broadcasting thread that it can continue.
+     * But if we continued due to cond_signal, we do not have to do that
+     * because the signaling thread knows that only one waiter continued.
+     */
+    EnterCriticalSection(&cond->waiters_lock);
+    cond->waiters--;
+    last_waiter = cond->was_broadcast && cond->waiters == 0;
+    LeaveCriticalSection(&cond->waiters_lock);
+
+    if (last_waiter) {
+        /*
+         * cond_broadcast was issued while mutex was held. This means
+         * that all other waiters have continued, but are contending
+         * for the mutex at the end of this function because the
+         * broadcasting thread did not leave cond_broadcast, yet.
+         * (This is so that it can be sure that each waiter has
+         * consumed exactly one slice of the semaphor.)
+         * The last waiter must tell the broadcasting thread that it
+         * can go on.
+         */
+        SetEvent(cond->continue_broadcast);
+        /*
+         * Now we go on to contend with all other waiters for
+         * the mutex. Auf in den Kampf!
+         */
+    }
+    /* lock external mutex again */
+    EnterCriticalSection(mutex);
+
+    return 0;
+}
+
+/*
+ * IMPORTANT: This implementation requires that pthread_cond_signal
+ * is called while the mutex is held that is used in the corresponding
+ * pthread_cond_wait calls!
+ */
+int pthread_cond_signal(pthread_cond_t *cond)
+{
+    int have_waiters;
+
+    EnterCriticalSection(&cond->waiters_lock);
+    have_waiters = cond->waiters > 0;
+    LeaveCriticalSection(&cond->waiters_lock);
+
+    /*
+     * Signal only when there are waiters
+     */
+    if (have_waiters)
+        return ReleaseSemaphore(cond->sema, 1, NULL) ?
+            0 : GetLastError();
+    else
+        return 0;
+}
+
+/*
+ * DOUBLY IMPORTANT: This implementation requires that pthread_cond_broadcast
+ * is called while the mutex is held that is used in the corresponding
+ * pthread_cond_wait calls!
+ */
+int pthread_cond_broadcast(pthread_cond_t *cond)
+{
+    EnterCriticalSection(&cond->waiters_lock);
+
+    if ((cond->was_broadcast = cond->waiters > 0)) {
+        /* wake up all waiters */
+        ReleaseSemaphore(cond->sema, cond->waiters, NULL);
+        LeaveCriticalSection(&cond->waiters_lock);
+        /*
+         * At this point all waiters continue. Each one takes its
+         * slice of the semaphor. Now it's our turn to wait: Since
+         * the external mutex is held, no thread can leave cond_wait,
+         * yet. For this reason, we can be sure that no thread gets
+         * a chance to eat *more* than one slice. OTOH, it means
+         * that the last waiter must send us a wake-up.
+         */
+        WaitForSingleObject(cond->continue_broadcast, INFINITE);
+        /*
+         * Since the external mutex is held, no thread can enter
+         * cond_wait, and, hence, it is safe to reset this flag
+         * without cond->waiters_lock held.
+         */
+        cond->was_broadcast = 0;
+    } else {
+        LeaveCriticalSection(&cond->waiters_lock);
+    }
+    return 0;
+}

numexpr/win32/pthread.h

@@ -0,0 +1,119 @@
+/*
+ * Code for simulating pthreads API on Windows.  This is Git-specific,
+ * but it is enough for Numexpr needs too.
+ *
+ * Copyright (C) 2009 Andrzej K. Haczewski <ahaczewski@gmail.com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ *
+ * DISCLAIMER: The implementation is Git-specific, it is subset of original
+ * Pthreads API, without lots of other features that Git doesn't use.
+ * Git also makes sure that the passed arguments are valid, so there's
+ * no need for double-checking.
+ */
+
+#ifndef PTHREAD_H
+#define PTHREAD_H
+
+#ifndef WIN32_LEAN_AND_MEAN
+#define WIN32_LEAN_AND_MEAN
+#endif
+
+#include <windows.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/*
+ * Defines that adapt Windows API threads to pthreads API
+ */
+#define pthread_mutex_t CRITICAL_SECTION
+
+#define pthread_mutex_init(a,b) InitializeCriticalSection((a))
+#define pthread_mutex_destroy(a) DeleteCriticalSection((a))
+#define pthread_mutex_lock EnterCriticalSection
+#define pthread_mutex_unlock LeaveCriticalSection
+
+/*
+ * Implement simple condition variable for Windows threads, based on ACE
+ * implementation.
+ *
+ * See original implementation: http://bit.ly/1vkDjo
+ * ACE homepage: http://www.cse.wustl.edu/~schmidt/ACE.html
+ * See also: http://www.cse.wustl.edu/~schmidt/win32-cv-1.html
+ */
+typedef struct {
+	LONG waiters;
+	int was_broadcast;
+	CRITICAL_SECTION waiters_lock;
+	HANDLE sema;
+	HANDLE continue_broadcast;
+} pthread_cond_t;
+
+extern int pthread_cond_init(pthread_cond_t *cond, const void *unused);
+extern int pthread_cond_destroy(pthread_cond_t *cond);
+extern int pthread_cond_wait(pthread_cond_t *cond, CRITICAL_SECTION *mutex);
+extern int pthread_cond_signal(pthread_cond_t *cond);
+extern int pthread_cond_broadcast(pthread_cond_t *cond);
+
+/*
+ * Simple thread creation implementation using pthread API
+ */
+typedef struct {
+	HANDLE handle;
+	void *(*start_routine)(void*);
+	void *arg;
+} pthread_t;
+
+extern int pthread_create(pthread_t *thread, const void *unused,
+			  void *(*start_routine)(void*), void *arg);
+
+/*
+ * To avoid the need of copying a struct, we use small macro wrapper to pass
+ * pointer to win32_pthread_join instead.
+ */
+#define pthread_join(a, b) win32_pthread_join(&(a), (b))
+
+extern int win32_pthread_join(pthread_t *thread, void **value_ptr);
+
+/*
+ * The POSIX signal system has a more developed interface than what's in
+ * Windows. We create a no-op shim layer to proivde enough of the API to
+ * pretend to support what's used when creating threads on POSIX systems.
+ */
+typedef int sigset_t;
+enum sigop {
+    SIG_BLOCK,
+    SIG_UNBLOCK,
+    SIG_SETMASK
+};
+
+static inline int sigemptyset(sigset_t *sigs) { return 0; }
+static inline int sigfillset(sigset_t *sigs) { return 0; }
+static inline int sigaddset(sigset_t *sigs, int sig) { return 0; }
+static inline int sigdelset(sigset_t *sigs, int sig) { return 0; }
+static inline int pthread_sigmask(int how, sigset_t *newmask,
+                                  sigset_t *oldmask) { return 0; }
+
+#ifdef __cplusplus
+} // extern "C"
+#endif
+
+#endif /* PTHREAD_H */

numexpr/win32/stdint.h

@@ -0,0 +1,235 @@
+/* ISO C9x  7.18  Integer types <stdint.h>
+ * Based on ISO/IEC SC22/WG14 9899 Committee draft (SC22 N2794)
+ *
+ *  THIS SOFTWARE IS NOT COPYRIGHTED
+ *
+ *  Contributor: Danny Smith <danny_r_smith_2001@yahoo.co.nz>
+ *
+ *  This source code is offered for use in the public domain. You may
+ *  use, modify or distribute it freely.
+ *
+ *  This code is distributed in the hope that it will be useful but
+ *  WITHOUT ANY WARRANTY. ALL WARRANTIES, EXPRESS OR IMPLIED ARE HEREBY
+ *  DISCLAIMED. This includes but is not limited to warranties of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *  Date: 2000-12-02
+ *
+ * mwb: This was modified in the following ways:
+ *
+ *      - make it compatible with Visual C++ 6 (which uses
+ *          non-standard keywords and suffixes for 64-bit types)
+ *      - some environments need stddef.h included (for wchar stuff?)
+ *      - handle the fact that Microsoft's limits.h header defines
+ *          SIZE_MAX
+ *      - make corrections for SIZE_MAX, INTPTR_MIN, INTPTR_MAX, UINTPTR_MAX,
+ *          PTRDIFF_MIN, PTRDIFF_MAX, SIG_ATOMIC_MIN, and SIG_ATOMIC_MAX
+ *          to be 64-bit aware.
+ */
+
+
+#ifndef _STDINT_H
+#define _STDINT_H
+#define __need_wint_t
+#define __need_wchar_t
+#include <wchar.h>
+#include <stddef.h>
+
+#if _MSC_VER && (_MSC_VER < 1300)
+/* using MSVC 6 or earlier - no "long long" type, but might have _int64 type */
+#define __STDINT_LONGLONG           __int64
+#define __STDINT_LONGLONG_SUFFIX    i64
+#else
+#define __STDINT_LONGLONG           long long
+#define __STDINT_LONGLONG_SUFFIX    LL
+#endif
+
+#if !defined( PASTE)
+#define PASTE2( x, y) x##y
+#define PASTE( x, y)  PASTE2( x, y)
+#endif /* PASTE */
+
+
+/* 7.18.1.1  Exact-width integer types */
+typedef signed char int8_t;
+typedef unsigned char   uint8_t;
+typedef short  int16_t;
+typedef unsigned short  uint16_t;
+typedef int  int32_t;
+typedef unsigned   uint32_t;
+typedef __STDINT_LONGLONG  int64_t;
+typedef unsigned __STDINT_LONGLONG   uint64_t;
+
+/* 7.18.1.2  Minimum-width integer types */
+typedef signed char int_least8_t;
+typedef unsigned char   uint_least8_t;
+typedef short  int_least16_t;
+typedef unsigned short  uint_least16_t;
+typedef int  int_least32_t;
+typedef unsigned   uint_least32_t;
+typedef __STDINT_LONGLONG  int_least64_t;
+typedef unsigned __STDINT_LONGLONG   uint_least64_t;
+
+/*  7.18.1.3  Fastest minimum-width integer types
+ *  Not actually guaranteed to be fastest for all purposes
+ *  Here we use the exact-width types for 8 and 16-bit ints.
+ */
+typedef char int_fast8_t;
+typedef unsigned char uint_fast8_t;
+typedef short  int_fast16_t;
+typedef unsigned short  uint_fast16_t;
+typedef int  int_fast32_t;
+typedef unsigned  int  uint_fast32_t;
+typedef __STDINT_LONGLONG  int_fast64_t;
+typedef unsigned __STDINT_LONGLONG   uint_fast64_t;
+
+/* 7.18.1.4  Integer types capable of holding object pointers */
+#ifndef _INTPTR_T_DEFINED
+#define _INTPTR_T_DEFINED
+#ifdef _WIN64
+typedef __STDINT_LONGLONG intptr_t
+#else
+typedef int intptr_t;
+#endif /* _WIN64 */
+#endif /* _INTPTR_T_DEFINED */
+
+#ifndef _UINTPTR_T_DEFINED
+#define _UINTPTR_T_DEFINED
+#ifdef _WIN64
+typedef unsigned __STDINT_LONGLONG uintptr_t
+#else
+typedef unsigned int uintptr_t;
+#endif /* _WIN64 */
+#endif /* _UINTPTR_T_DEFINED */
+
+/* 7.18.1.5  Greatest-width integer types */
+typedef __STDINT_LONGLONG  intmax_t;
+typedef unsigned __STDINT_LONGLONG   uintmax_t;
+
+/* 7.18.2  Limits of specified-width integer types */
+#if !defined ( __cplusplus) || defined (__STDC_LIMIT_MACROS)
+
+/* 7.18.2.1  Limits of exact-width integer types */
+#define INT8_MIN (-128)
+#define INT16_MIN (-32768)
+#define INT32_MIN (-2147483647 - 1)
+#define INT64_MIN  (PASTE( -9223372036854775807, __STDINT_LONGLONG_SUFFIX) - 1)
+
+#define INT8_MAX 127
+#define INT16_MAX 32767
+#define INT32_MAX 2147483647
+#define INT64_MAX (PASTE( 9223372036854775807, __STDINT_LONGLONG_SUFFIX))
+
+#define UINT8_MAX 0xff /* 255U */
+#define UINT16_MAX 0xffff /* 65535U */
+#define UINT32_MAX 0xffffffff  /* 4294967295U */
+#define UINT64_MAX (PASTE( 0xffffffffffffffffU, __STDINT_LONGLONG_SUFFIX)) /* 18446744073709551615ULL */
+
+/* 7.18.2.2  Limits of minimum-width integer types */
+#define INT_LEAST8_MIN INT8_MIN
+#define INT_LEAST16_MIN INT16_MIN
+#define INT_LEAST32_MIN INT32_MIN
+#define INT_LEAST64_MIN INT64_MIN
+
+#define INT_LEAST8_MAX INT8_MAX
+#define INT_LEAST16_MAX INT16_MAX
+#define INT_LEAST32_MAX INT32_MAX
+#define INT_LEAST64_MAX INT64_MAX
+
+#define UINT_LEAST8_MAX UINT8_MAX
+#define UINT_LEAST16_MAX UINT16_MAX
+#define UINT_LEAST32_MAX UINT32_MAX
+#define UINT_LEAST64_MAX UINT64_MAX
+
+/* 7.18.2.3  Limits of fastest minimum-width integer types */
+#define INT_FAST8_MIN INT8_MIN
+#define INT_FAST16_MIN INT16_MIN
+#define INT_FAST32_MIN INT32_MIN
+#define INT_FAST64_MIN INT64_MIN
+
+#define INT_FAST8_MAX INT8_MAX
+#define INT_FAST16_MAX INT16_MAX
+#define INT_FAST32_MAX INT32_MAX
+#define INT_FAST64_MAX INT64_MAX
+
+#define UINT_FAST8_MAX UINT8_MAX
+#define UINT_FAST16_MAX UINT16_MAX
+#define UINT_FAST32_MAX UINT32_MAX
+#define UINT_FAST64_MAX UINT64_MAX
+
+/* 7.18.2.4  Limits of integer types capable of holding
+    object pointers */
+#ifdef _WIN64
+#define INTPTR_MIN INT64_MIN
+#define INTPTR_MAX INT64_MAX
+#define UINTPTR_MAX UINT64_MAX
+#else
+#define INTPTR_MIN INT32_MIN
+#define INTPTR_MAX INT32_MAX
+#define UINTPTR_MAX UINT32_MAX
+#endif /* _WIN64 */
+
+/* 7.18.2.5  Limits of greatest-width integer types */
+#define INTMAX_MIN INT64_MIN
+#define INTMAX_MAX INT64_MAX
+#define UINTMAX_MAX UINT64_MAX
+
+/* 7.18.3  Limits of other integer types */
+#define PTRDIFF_MIN INTPTR_MIN
+#define PTRDIFF_MAX INTPTR_MAX
+
+#define SIG_ATOMIC_MIN INTPTR_MIN
+#define SIG_ATOMIC_MAX INTPTR_MAX
+
+/* we need to check for SIZE_MAX already defined because MS defines it in limits.h */
+#ifndef SIZE_MAX
+#define SIZE_MAX UINTPTR_MAX
+#endif
+
+#ifndef WCHAR_MIN  /* also in wchar.h */
+#define WCHAR_MIN 0
+#define WCHAR_MAX ((wchar_t)-1) /* UINT16_MAX */
+#endif
+
+/*
+ * wint_t is unsigned short for compatibility with MS runtime
+ */
+#define WINT_MIN 0
+#define WINT_MAX ((wint_t)-1) /* UINT16_MAX */
+
+#endif /* !defined ( __cplusplus) || defined __STDC_LIMIT_MACROS */
+
+
+/* 7.18.4  Macros for integer constants */
+#if !defined ( __cplusplus) || defined (__STDC_CONSTANT_MACROS)
+
+/* 7.18.4.1  Macros for minimum-width integer constants
+
+    Accoding to Douglas Gwyn <gwyn@arl.mil>:
+	"This spec was changed in ISO/IEC 9899:1999 TC1; in ISO/IEC
+	9899:1999 as initially published, the expansion was required
+	to be an integer constant of precisely matching type, which
+	is impossible to accomplish for the shorter types on most
+	platforms, because C99 provides no standard way to designate
+	an integer constant with width less than that of type int.
+	TC1 changed this to require just an integer constant
+	*expression* with *promoted* type."
+*/
+
+#define INT8_C(val) ((int8_t) + (val))
+#define UINT8_C(val) ((uint8_t) + (val##U))
+#define INT16_C(val) ((int16_t) + (val))
+#define UINT16_C(val) ((uint16_t) + (val##U))
+
+#define INT32_C(val) val##L
+#define UINT32_C(val) val##UL
+#define INT64_C(val) (PASTE( val, __STDINT_LONGLONG_SUFFIX))
+#define UINT64_C(val)(PASTE( PASTE( val, U), __STDINT_LONGLONG_SUFFIX))
+
+/* 7.18.4.2  Macros for greatest-width integer constants */
+#define INTMAX_C(val)  INT64_C(val)
+#define UINTMAX_C(val) UINT64_C(val)
+
+#endif  /* !defined ( __cplusplus) || defined __STDC_CONSTANT_MACROS */
+
+#endif