make optional C99 double routines mandatory

3 files changed, 13 insertions, 368 deletions

diff --git a/numpy/core/setup.py b/numpy/core/setup.py
index 9e05aded0..571d067e0 100644
--- a/numpy/core/setup.py
+++ b/numpy/core/setup.py
@@ -145,7 +145,8 @@ def check_math_capabilities(config, ext, moredefs, mathlibs):
             headers=headers,
         )
 
-    def check_funcs_once(funcs_name, headers=["feature_detection_math.h"]):
+    def check_funcs_once(funcs_name, headers=["feature_detection_math.h"],
+                         add_to_moredefs=True):
         call = dict([(f, True) for f in funcs_name])
         call_args = dict([(f, FUNC_CALL_ARGS[f]) for f in funcs_name])
         st = config.check_funcs_once(
@@ -156,7 +157,7 @@ def check_math_capabilities(config, ext, moredefs, mathlibs):
             call_args=call_args,
             headers=headers,
         )
-        if st:
+        if st and add_to_moredefs:
             moredefs.extend([(fname2def(f), 1) for f in funcs_name])
         return st
 
@@ -173,7 +174,7 @@ def check_math_capabilities(config, ext, moredefs, mathlibs):
             return 1
 
     #use_msvc = config.check_decl("_MSC_VER")
-    if not check_funcs_once(MANDATORY_FUNCS):
+    if not check_funcs_once(MANDATORY_FUNCS, add_to_moredefs=False):
         raise SystemError("One of the required function to build numpy is not"
                 " available (the list is %s)." % str(MANDATORY_FUNCS))
 
@@ -184,20 +185,12 @@ def check_math_capabilities(config, ext, moredefs, mathlibs):
     # config.h in the public namespace, so we have a clash for the common
     # functions we test. We remove every function tested by python's
     # autoconf, hoping their own test are correct
-    for f in OPTIONAL_STDFUNCS_MAYBE:
-        if config.check_decl(fname2def(f),
-                    headers=["Python.h", "math.h"]):
-            if f in OPTIONAL_STDFUNCS:
-                OPTIONAL_STDFUNCS.remove(f)
-            else:
-                OPTIONAL_FILE_FUNCS.remove(f)
+    for f in OPTIONAL_FUNCS_MAYBE:
+        if config.check_decl(fname2def(f), headers=["Python.h"]):
+            OPTIONAL_FILE_FUNCS.remove(f)
 
-
-    check_funcs(OPTIONAL_STDFUNCS)
     check_funcs(OPTIONAL_FILE_FUNCS, headers=["feature_detection_stdio.h"])
     check_funcs(OPTIONAL_MISC_FUNCS, headers=["feature_detection_misc.h"])
-    
-
 
     for h in OPTIONAL_HEADERS:
         if config.check_func("", decl=False, call=False, headers=[h]):
diff --git a/numpy/core/setup_common.py b/numpy/core/setup_common.py
index a751acc6f..d618fac69 100644
--- a/numpy/core/setup_common.py
+++ b/numpy/core/setup_common.py
@@ -124,9 +124,8 @@ MANDATORY_FUNCS = ["sin", "cos", "tan", "sinh", "cosh", "tanh", "fabs",
         "floor", "ceil", "sqrt", "log10", "log", "exp", "asin",
         "acos", "atan", "fmod", 'modf', 'frexp', 'ldexp']
 
-# Standard functions which may not be available and for which we have a
-# replacement implementation. Note that some of these are C99 functions.
-OPTIONAL_STDFUNCS = ["expm1", "log1p", "acosh", "asinh", "atanh",
+# Some of these are C99 functions.
+MANDATORY_FUNCS += ["expm1", "log1p", "acosh", "asinh", "atanh",
         "rint", "trunc", "exp2", "log2", "hypot", "atan2", "pow",
         "copysign", "nextafter", "strtoll", "strtoull", "cbrt"]
 
@@ -227,9 +226,8 @@ OPTIONAL_FUNCTION_ATTRIBUTES_WITH_INTRINSICS = [('__attribute__((target("avx2,fm
 # variable attributes tested via "int %s a" % attribute
 OPTIONAL_VARIABLE_ATTRIBUTES = ["__thread", "__declspec(thread)"]
 
-# Subset of OPTIONAL_STDFUNCS which may already have HAVE_* defined by Python.h
-OPTIONAL_STDFUNCS_MAYBE = [
-    "expm1", "log1p", "acosh", "atanh", "asinh", "hypot", "copysign",
+# Subset of OPTIONAL_*_FUNCS which may already have HAVE_* defined by Python.h
+OPTIONAL_FUNCS_MAYBE = [
     "ftello", "fseeko"
     ]
 
diff --git a/numpy/core/src/npymath/npy_math_internal.h.src b/numpy/core/src/npymath/npy_math_internal.h.src
index 0364aabb1..b5e415c8a 100644
--- a/numpy/core/src/npymath/npy_math_internal.h.src
+++ b/numpy/core/src/npymath/npy_math_internal.h.src
@@ -88,297 +88,13 @@ static const npy_uint64 MAGIC64[] = {0x5555555555555555ull, 0x3333333333333333ul
  */
 
 /* Original code by Konrad Hinsen.  */
-#ifndef HAVE_EXPM1
-NPY_INPLACE double npy_expm1(double x)
-{
-    if (npy_isinf(x) && x > 0) {
-        return x;
-    }
-    else {
-        const double u = npy_exp(x);
-
-        if (u == 1.0) {
-            return x;
-        } else if (u - 1.0 == -1.0) {
-            return -1;
-        } else {
-            return (u - 1.0) * x/npy_log(u);
-        }
-    }
-}
-#endif
-
-#ifndef HAVE_LOG1P
-NPY_INPLACE double npy_log1p(double x)
-{
-    if (npy_isinf(x) && x > 0) {
-        return x;
-    }
-    else {
-        const double u = 1. + x;
-        const double d = u - 1.;
-
-        if (d == 0) {
-            return x;
-        } else {
-            return npy_log(u) * x / d;
-        }
-    }
-}
-#endif
-
 /* Taken from FreeBSD mlib, adapted for numpy
  *
  * XXX: we could be a bit faster by reusing high/low words for inf/nan
  * classification instead of calling npy_isinf/npy_isnan: we should have some
  * macros for this, though, instead of doing it manually
  */
-#ifndef HAVE_ATAN2
-/* XXX: we should have this in npy_math.h */
-#define NPY_DBL_EPSILON 1.2246467991473531772E-16
-NPY_INPLACE double npy_atan2(double y, double x)
-{
-    npy_int32 k, m, iy, ix, hx, hy;
-    npy_uint32 lx,ly;
-    double z;
-
-    EXTRACT_WORDS(hx, lx, x);
-    ix = hx & 0x7fffffff;
-    EXTRACT_WORDS(hy, ly, y);
-    iy = hy & 0x7fffffff;
-
-    /* if x or y is nan, return nan */
-    if (npy_isnan(x * y)) {
-        return x + y;
-    }
-
-    if (x == 1.0) {
-        return npy_atan(y);
-    }
-
-    m = 2 * (npy_signbit((x)) != 0) + (npy_signbit((y)) != 0);
-    if (y == 0.0) {
-        switch(m) {
-        case 0:
-        case 1: return  y;  /* atan(+-0,+anything)=+-0 */
-        case 2: return  NPY_PI;/* atan(+0,-anything) = pi */
-        case 3: return -NPY_PI;/* atan(-0,-anything) =-pi */
-        }
-    }
-
-    if (x == 0.0) {
-        return y > 0 ? NPY_PI_2 : -NPY_PI_2;
-    }
-
-    if (npy_isinf(x)) {
-        if (npy_isinf(y)) {
-            switch(m) {
-                case 0: return  NPY_PI_4;/* atan(+INF,+INF) */
-                case 1: return -NPY_PI_4;/* atan(-INF,+INF) */
-                case 2: return  3.0*NPY_PI_4;/*atan(+INF,-INF)*/
-                case 3: return -3.0*NPY_PI_4;/*atan(-INF,-INF)*/
-            }
-        } else {
-            switch(m) {
-                case 0: return  NPY_PZERO;  /* atan(+...,+INF) */
-                case 1: return  NPY_NZERO;  /* atan(-...,+INF) */
-                case 2: return  NPY_PI;  /* atan(+...,-INF) */
-                case 3: return -NPY_PI;  /* atan(-...,-INF) */
-            }
-        }
-    }
-
-    if (npy_isinf(y)) {
-        return y > 0 ? NPY_PI_2 : -NPY_PI_2;
-    }
-
-    /* compute y/x */
-    k = (iy - ix) >> 20;
-    if (k > 60) {            /* |y/x| >  2**60 */
-        z = NPY_PI_2 + 0.5 * NPY_DBL_EPSILON;
-        m &= 1;
-    } else if (hx < 0 && k < -60) {
-        z = 0.0;    /* 0 > |y|/x > -2**-60 */
-    } else {
-        z = npy_atan(npy_fabs(y/x));        /* safe to do y/x */
-    }
-
-    switch (m) {
-        case 0: return  z  ;    /* atan(+,+) */
-        case 1: return -z  ;    /* atan(-,+) */
-        case 2: return  NPY_PI - (z - NPY_DBL_EPSILON);/* atan(+,-) */
-        default: /* case 3 */
-            return  (z - NPY_DBL_EPSILON) - NPY_PI;/* atan(-,-) */
-    }
-}
-
-#endif
-
-#ifndef HAVE_HYPOT
-NPY_INPLACE double npy_hypot(double x, double y)
-{
-    double yx;
-
-    if (npy_isinf(x) || npy_isinf(y)) {
-        return NPY_INFINITY;
-    }
-
-    if (npy_isnan(x) || npy_isnan(y)) {
-        return NPY_NAN;
-    }
-
-    x = npy_fabs(x);
-    y = npy_fabs(y);
-    if (x < y) {
-        double temp = x;
-        x = y;
-        y = temp;
-    }
-    if (x == 0.) {
-        return 0.;
-    }
-    else {
-        yx = y/x;
-        return x*npy_sqrt(1.+yx*yx);
-    }
-}
-#endif
-
-#ifndef HAVE_ACOSH
-NPY_INPLACE double npy_acosh(double x)
-{
-    if (x < 1.0) {
-        return NPY_NAN;
-    }
-
-    if (npy_isfinite(x)) {
-        if (x > 1e8) {
-             return npy_log(x) + NPY_LOGE2;
-        }
-        else {
-            double u = x - 1.0;
-            return npy_log1p(u + npy_sqrt(2*u + u*u));
-        }
-    }
-    return x;
-}
-#endif
-
-#ifndef HAVE_ASINH
-NPY_INPLACE double npy_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 = npy_sqrt(x*x + 1);
-    }
-    return sign*npy_log1p(x*(1.0 + x/(d+1)));
-}
-#endif
-
-#ifndef HAVE_ATANH
-NPY_INPLACE double npy_atanh(double x)
-{
-    if (x > 0) {
-        return -0.5*npy_log1p(-2.0*x/(1.0 + x));
-    }
-    else {
-        return 0.5*npy_log1p(2.0*x/(1.0 - x));
-    }
-}
-#endif
-
-#ifndef HAVE_RINT
-#if defined(_MSC_VER) && (_MSC_VER == 1500) && !defined(_WIN64)
-#pragma optimize("", off)
-#endif
-NPY_INPLACE double npy_rint(double x)
-{
-    double y, r;
-
-    y = npy_floor(x);
-    r = x - y;
-
-    if (r > 0.5) {
-        y += 1.0;
-    }
-
-    /* Round to nearest even */
-    if (r == 0.5) {
-        r = y - 2.0*npy_floor(0.5*y);
-        if (r == 1.0) {
-            y += 1.0;
-        }
-    }
-    return y;
-}
-#if defined(_MSC_VER) && (_MSC_VER == 1500) && !defined(_WIN64)
-#pragma optimize("", on)
-#endif
-#endif
-
-#ifndef HAVE_TRUNC
-NPY_INPLACE double npy_trunc(double x)
-{
-    return x < 0 ? npy_ceil(x) : npy_floor(x);
-}
-#endif
-
-#ifndef HAVE_EXP2
-NPY_INPLACE double npy_exp2(double x)
-{
-    return npy_exp(NPY_LOGE2*x);
-}
-#endif
-
-#ifndef HAVE_LOG2
-NPY_INPLACE double npy_log2(double x)
-{
-#ifdef HAVE_FREXP
-    if (!npy_isfinite(x) || x <= 0.) {
-        /* special value result */
-        return npy_log(x);
-    }
-    else {
-        /*
-         * fallback implementation copied from python3.4 math.log2
-         * provides int(log(2**i)) == i for i 1-64 in default rounding mode.
-         *
-         * We want log2(m * 2**e) == log(m) / log(2) + e.  Care is needed when
-         * x is just greater than 1.0: in that case e is 1, log(m) is negative,
-         * and we get significant cancellation error from the addition of
-         * log(m) / log(2) to e.  The slight rewrite of the expression below
-         * avoids this problem.
-         */
-        int e;
-        double m = frexp(x, &e);
-        if (x >= 1.0) {
-            return log(2.0 * m) / log(2.0) + (e - 1);
-        }
-        else {
-            return log(m) / log(2.0) + e;
-        }
-    }
-#else
-    /* does not provide int(log(2**i)) == i */
-    return NPY_LOG2E * npy_log(x);
-#endif
-}
-#endif
-
 /*
- * if C99 extensions not available then define dummy functions that use the
- * double versions for
  *
  * sin, cos, tan
  * sinh, cosh, tanh,
@@ -388,7 +104,7 @@ NPY_INPLACE double npy_log2(double x)
  * asinh, acosh, atanh
  *
  * hypot, atan2, pow, fmod, modf
- * ldexp, frexp
+ * ldexp, frexp, cbrt
  *
  * We assume the above are always available in their double versions.
  *
@@ -404,40 +120,6 @@ NPY_INPLACE double npy_log2(double x)
  * #C = L,F#
  */
 
-/**begin repeat1
- * #kind = sin,cos,tan,sinh,cosh,tanh,fabs,floor,ceil,rint,trunc,sqrt,log10,
- *         log,exp,expm1,asin,acos,atan,asinh,acosh,atanh,log1p,exp2,log2#
- * #KIND = SIN,COS,TAN,SINH,COSH,TANH,FABS,FLOOR,CEIL,RINT,TRUNC,SQRT,LOG10,
- *         LOG,EXP,EXPM1,ASIN,ACOS,ATAN,ASINH,ACOSH,ATANH,LOG1P,EXP2,LOG2#
- */
-
-#ifdef @kind@@c@
-#undef @kind@@c@
-#endif
-#ifndef HAVE_@KIND@@C@
-NPY_INPLACE @type@ npy_@kind@@c@(@type@ x)
-{
-    return (@type@) npy_@kind@((double)x);
-}
-#endif
-
-/**end repeat1**/
-
-/**begin repeat1
- * #kind = atan2,hypot,pow,fmod,copysign#
- * #KIND = ATAN2,HYPOT,POW,FMOD,COPYSIGN#
- */
-#ifdef @kind@@c@
-#undef @kind@@c@
-#endif
-#ifndef HAVE_@KIND@@C@
-NPY_INPLACE @type@ npy_@kind@@c@(@type@ x, @type@ y)
-{
-    return (@type@) npy_@kind@((double)x, (double) y);
-}
-#endif
-/**end repeat1**/
-
 #ifdef modf@c@
 #undef modf@c@
 #endif
@@ -503,11 +185,8 @@ NPY_INPLACE @type@ npy_frexp@c@(@type@ x, int* exp)
 /**begin repeat1
  * #kind = sin,cos,tan,sinh,cosh,tanh,fabs,floor,ceil,rint,trunc,sqrt,log10,
  *         log,exp,expm1,asin,acos,atan,asinh,acosh,atanh,log1p,exp2,log2#
- * #KIND = SIN,COS,TAN,SINH,COSH,TANH,FABS,FLOOR,CEIL,RINT,TRUNC,SQRT,LOG10,
- *         LOG,EXP,EXPM1,ASIN,ACOS,ATAN,ASINH,ACOSH,ATANH,LOG1P,EXP2,LOG2#
  * #TRIG_WORKAROUND = WORKAROUND_APPLE_TRIG_BUG*3, 0*22#
  */
-#ifdef HAVE_@KIND@@C@
 NPY_INPLACE @type@ npy_@kind@@c@(@type@ x)
 {
 #if @TRIG_WORKAROUND@
@@ -517,7 +196,6 @@ NPY_INPLACE @type@ npy_@kind@@c@(@type@ x)
 #endif
     return NPY__FP_SFX(@kind@)(x);
 }
-#endif
 
 /**end repeat1**/
 
@@ -527,57 +205,33 @@ NPY_INPLACE @type@ npy_@kind@@c@(@type@ x)
  * #kind = atan2,hypot,pow,fmod,copysign#
  * #KIND = ATAN2,HYPOT,POW,FMOD,COPYSIGN#
  */
-#ifdef HAVE_@KIND@@C@
 NPY_INPLACE @type@ npy_@kind@@c@(@type@ x, @type@ y)
 {
     return NPY__FP_SFX(@kind@)(x, y);
 }
-#endif
 /**end repeat1**/
 
-#ifdef HAVE_MODF@C@
 NPY_INPLACE @type@ npy_modf@c@(@type@ x, @type@ *iptr)
 {
     return NPY__FP_SFX(modf)(x, iptr);
 }
-#endif
 
-#ifdef HAVE_LDEXP@C@
 NPY_INPLACE @type@ npy_ldexp@c@(@type@ x, int exp)
 {
     return NPY__FP_SFX(ldexp)(x, exp);
 }
-#endif
 
-#ifdef HAVE_FREXP@C@
 NPY_INPLACE @type@ npy_frexp@c@(@type@ x, int* exp)
 {
     return NPY__FP_SFX(frexp)(x, exp);
 }
-#endif
 
-/* C99 but not mandatory */
+/* C99 mandatory */
 
-#ifndef HAVE_CBRT@C@
-NPY_INPLACE @type@ npy_cbrt@c@(@type@ x)
-{
-    /* don't set invalid flag */
-    if (npy_isnan(x)) {
-        return NPY_NAN;
-    }
-    else if (x < 0) {
-        return -npy_pow@c@(-x, 1. / 3.);
-    }
-    else {
-        return npy_pow@c@(x, 1. / 3.);
-    }
-}
-#else
 NPY_INPLACE @type@ npy_cbrt@c@(@type@ x)
 {
     return NPY__FP_SFX(cbrt)(x);
 }
-#endif
 #undef NPY__FP_SFX
 /**end repeat**/