MAINT: Translate ieee754.c.src to C++ using templates. (#21367)

* [Road to C++] Rename ieee754.c.src into ieee754.cpp.src

* [Road to C++] Remove all use of internal template engine

Replace them by regular C++ templates.

* [Road to C++] mv ieee754.cpp.src ro ieee754.cpp

Clang-format it, and update build system accordingly.

* Fix 'unnecessary parentheses in declaration of ‘__NPY_UNUSED_TAGGEDc’ [-Wparentheses]'

* [Road to C++] Fix misc conversion issues in ieee754.cpp

- hex float literal are not supported in C++
- avoid overload issues when numpy_longdouble ultimately is equal to double
- _statusfp2 exepcts unsigned arguments

1 files changed, 339 insertions, 237 deletions

diff --git a/numpy/core/src/npymath/ieee754.c.src b/numpy/core/src/npymath/ieee754.cpp
index 4e6ddb712..81573d97a 100644
--- a/numpy/core/src/npymath/ieee754.c.src
+++ b/numpy/core/src/npymath/ieee754.cpp
@@ -4,12 +4,14 @@
  *
  * Low-level routines related to IEEE-754 format
  */
+#include "numpy/utils.h"
+
 #include "npy_math_common.h"
 #include "npy_math_private.h"
-#include "numpy/utils.h"
 
 #ifndef HAVE_COPYSIGN
-double npy_copysign(double x, double y)
+double
+npy_copysign(double x, double y)
 {
     npy_uint32 hx, hy;
     GET_HIGH_WORD(hx, x);
@@ -20,7 +22,8 @@ double npy_copysign(double x, double y)
 #endif
 
 /*
- The below code is provided for compilers which do not yet provide C11 compatibility (gcc 4.5 and older)
+ The below code is provided for compilers which do not yet provide C11
+ compatibility (gcc 4.5 and older)
  */
 #ifndef LDBL_TRUE_MIN
 #define LDBL_TRUE_MIN __LDBL_DENORM_MIN__
@@ -29,14 +32,16 @@ double npy_copysign(double x, double y)
 #if !defined(HAVE_DECL_SIGNBIT)
 #include "_signbit.c"
 
-int _npy_signbit_f(float x)
+int
+_npy_signbit_f(float x)
 {
-    return _npy_signbit_d((double) x);
+    return _npy_signbit_d((double)x);
 }
 
-int _npy_signbit_ld(long double x)
+int
+_npy_signbit_ld(long double x)
 {
-    return _npy_signbit_d((double) x);
+    return _npy_signbit_d((double)x);
 }
 #endif
 
@@ -47,44 +52,47 @@ int _npy_signbit_ld(long double x)
  * p >= 0, returnx x + nulp
  * p < 0, returnx x - nulp
  */
-static double _next(double x, int p)
+static double
+_next(double x, int p)
 {
     volatile double t;
     npy_int32 hx, hy, ix;
     npy_uint32 lx;
 
     EXTRACT_WORDS(hx, lx, x);
-    ix = hx & 0x7fffffff;       /* |x| */
+    ix = hx & 0x7fffffff; /* |x| */
 
-    if (((ix >= 0x7ff00000) && ((ix - 0x7ff00000) | lx) != 0))        /* x is nan */
+    if (((ix >= 0x7ff00000) && ((ix - 0x7ff00000) | lx) != 0)) /* x is nan */
         return x;
-    if ((ix | lx) == 0) {       /* x == 0 */
+    if ((ix | lx) == 0) { /* x == 0 */
         if (p >= 0) {
-            INSERT_WORDS(x, 0x0, 1);    /* return +minsubnormal */
-        } else {
-            INSERT_WORDS(x, 0x80000000, 1);    /* return -minsubnormal */
+            INSERT_WORDS(x, 0x0, 1); /* return +minsubnormal */
+        }
+        else {
+            INSERT_WORDS(x, 0x80000000, 1); /* return -minsubnormal */
         }
         t = x * x;
         if (t == x)
             return t;
         else
-            return x;           /* raise underflow flag */
+            return x; /* raise underflow flag */
     }
-    if (p < 0) {     /* x -= ulp */
+    if (p < 0) { /* x -= ulp */
         if (lx == 0)
             hx -= 1;
         lx -= 1;
-    } else {         /* x += ulp */
+    }
+    else { /* x += ulp */
         lx += 1;
         if (lx == 0)
             hx += 1;
     }
     hy = hx & 0x7ff00000;
     if (hy >= 0x7ff00000)
-        return x + x;           /* overflow  */
-    if (hy < 0x00100000) {      /* underflow */
+        return x + x;      /* overflow  */
+    if (hy < 0x00100000) { /* underflow */
         t = x * x;
-        if (t != x) {           /* raise underflow flag */
+        if (t != x) { /* raise underflow flag */
             INSERT_WORDS(x, hx, lx);
             return x;
         }
@@ -93,39 +101,42 @@ static double _next(double x, int p)
     return x;
 }
 
-static float _nextf(float x, int p)
+static float
+_next(float x, int p)
 {
     volatile float t;
     npy_int32 hx, hy, ix;
 
     GET_FLOAT_WORD(hx, x);
-    ix = hx & 0x7fffffff;       /* |x| */
+    ix = hx & 0x7fffffff; /* |x| */
 
-    if ((ix > 0x7f800000))      /* x is nan */
+    if ((ix > 0x7f800000)) /* x is nan */
         return x;
-    if (ix == 0) {              /* x == 0 */
+    if (ix == 0) { /* x == 0 */
         if (p >= 0) {
             SET_FLOAT_WORD(x, 0x0 | 1); /* return +minsubnormal */
-        } else {
+        }
+        else {
             SET_FLOAT_WORD(x, 0x80000000 | 1); /* return -minsubnormal */
         }
         t = x * x;
         if (t == x)
             return t;
         else
-            return x;           /* raise underflow flag */
+            return x; /* raise underflow flag */
     }
-    if (p < 0) {            /* x -= ulp */
+    if (p < 0) { /* x -= ulp */
         hx -= 1;
-    } else {                /* x += ulp */
+    }
+    else { /* x += ulp */
         hx += 1;
     }
     hy = hx & 0x7f800000;
     if (hy >= 0x7f800000)
-        return x + x;           /* overflow  */
-    if (hy < 0x00800000) {      /* underflow */
+        return x + x;      /* overflow  */
+    if (hy < 0x00800000) { /* underflow */
         t = x * x;
-        if (t != x) {           /* raise underflow flag */
+        if (t != x) { /* raise underflow flag */
             SET_FLOAT_WORD(x, hx);
             return x;
         }
@@ -135,129 +146,140 @@ static float _nextf(float x, int p)
 }
 
 #if defined(HAVE_LDOUBLE_DOUBLE_DOUBLE_BE) || \
-    defined(HAVE_LDOUBLE_DOUBLE_DOUBLE_LE)
+        defined(HAVE_LDOUBLE_DOUBLE_DOUBLE_LE)
 
 /*
  * FIXME: this is ugly and untested. The asm part only works with gcc, and we
  * should consolidate the GET_LDOUBLE* / SET_LDOUBLE macros
  */
-#define math_opt_barrier(x) \
-        ({ __typeof (x) __x = x; __asm ("" : "+m" (__x)); __x; })
-#define math_force_eval(x) __asm __volatile ("" : : "m" (x))
+#define math_opt_barrier(x)    \
+    ({                         \
+        __typeof(x) __x = x;   \
+        __asm("" : "+m"(__x)); \
+        __x;                   \
+    })
+#define math_force_eval(x) __asm __volatile("" : : "m"(x))
 
 /* only works for big endian */
-typedef union
-{
+typedef union {
     npy_longdouble value;
-    struct
-    {
+    struct {
         npy_uint64 msw;
         npy_uint64 lsw;
     } parts64;
-    struct
-    {
+    struct {
         npy_uint32 w0, w1, w2, w3;
     } parts32;
 } ieee854_long_double_shape_type;
 
 /* Get two 64 bit ints from a long double.  */
 
-#define GET_LDOUBLE_WORDS64(ix0,ix1,d) \
-do {                                   \
-  ieee854_long_double_shape_type qw_u; \
-  qw_u.value = (d);                    \
-  (ix0) = qw_u.parts64.msw;            \
-  (ix1) = qw_u.parts64.lsw;            \
-} while (0)
+#define GET_LDOUBLE_WORDS64(ix0, ix1, d)     \
+    do {                                     \
+        ieee854_long_double_shape_type qw_u; \
+        qw_u.value = (d);                    \
+        (ix0) = qw_u.parts64.msw;            \
+        (ix1) = qw_u.parts64.lsw;            \
+    } while (0)
 
 /* Set a long double from two 64 bit ints.  */
 
-#define SET_LDOUBLE_WORDS64(d,ix0,ix1) \
-do {                                   \
-  ieee854_long_double_shape_type qw_u; \
-  qw_u.parts64.msw = (ix0);            \
-  qw_u.parts64.lsw = (ix1);            \
-  (d) = qw_u.value;                    \
-} while (0)
+#define SET_LDOUBLE_WORDS64(d, ix0, ix1)     \
+    do {                                     \
+        ieee854_long_double_shape_type qw_u; \
+        qw_u.parts64.msw = (ix0);            \
+        qw_u.parts64.lsw = (ix1);            \
+        (d) = qw_u.value;                    \
+    } while (0)
 
-static npy_longdouble _nextl(npy_longdouble x, int p)
+static long double
+_next(long double x, int p)
 {
-    npy_int64 hx,ihx,ilx;
+    npy_int64 hx, ihx, ilx;
     npy_uint64 lx;
     npy_longdouble u;
+    const npy_longdouble eps = exp2l(-105.); // 0x1.0000000000000p-105L
 
     GET_LDOUBLE_WORDS64(hx, lx, x);
-    ihx = hx & 0x7fffffffffffffffLL;      /* |hx| */
-    ilx = lx & 0x7fffffffffffffffLL;      /* |lx| */
+    ihx = hx & 0x7fffffffffffffffLL; /* |hx| */
+    ilx = lx & 0x7fffffffffffffffLL; /* |lx| */
 
-    if(((ihx & 0x7ff0000000000000LL)==0x7ff0000000000000LL)&&
-       ((ihx & 0x000fffffffffffffLL)!=0)) {
+    if (((ihx & 0x7ff0000000000000LL) == 0x7ff0000000000000LL) &&
+        ((ihx & 0x000fffffffffffffLL) != 0)) {
         return x; /* signal the nan */
     }
-    if(ihx == 0 && ilx == 0) {          /* x == 0 */
-        SET_LDOUBLE_WORDS64(x, p, 0ULL);/* return +-minsubnormal */
+    if (ihx == 0 && ilx == 0) {          /* x == 0 */
+        SET_LDOUBLE_WORDS64(x, p, 0ULL); /* return +-minsubnormal */
         u = x * x;
         if (u == x) {
             return u;
-        } else {
-            return x;           /* raise underflow flag */
+        }
+        else {
+            return x; /* raise underflow flag */
         }
     }
 
-    if(p < 0) { /* p < 0, x -= ulp */
-        if((hx==0xffefffffffffffffLL)&&(lx==0xfc8ffffffffffffeLL))
-            return x+x; /* overflow, return -inf */
+    if (p < 0) { /* p < 0, x -= ulp */
+        if ((hx == 0xffefffffffffffffLL) && (lx == 0xfc8ffffffffffffeLL))
+            return x + x; /* overflow, return -inf */
         if (hx >= 0x7ff0000000000000LL) {
-            SET_LDOUBLE_WORDS64(u,0x7fefffffffffffffLL,0x7c8ffffffffffffeLL);
+            SET_LDOUBLE_WORDS64(u, 0x7fefffffffffffffLL, 0x7c8ffffffffffffeLL);
             return u;
         }
-        if(ihx <= 0x0360000000000000LL) {  /* x <= LDBL_MIN */
-            u = math_opt_barrier (x);
+        if (ihx <= 0x0360000000000000LL) { /* x <= LDBL_MIN */
+            u = math_opt_barrier(x);
             x -= LDBL_TRUE_MIN;
-            if (ihx < 0x0360000000000000LL
-                    || (hx > 0 && (npy_int64) lx <= 0)
-                    || (hx < 0 && (npy_int64) lx > 1)) {
+            if (ihx < 0x0360000000000000LL || (hx > 0 && (npy_int64)lx <= 0) ||
+                (hx < 0 && (npy_int64)lx > 1)) {
                 u = u * u;
-                math_force_eval (u);        /* raise underflow flag */
+                math_force_eval(u); /* raise underflow flag */
             }
             return x;
         }
         if (ihx < 0x06a0000000000000LL) { /* ulp will denormal */
-            SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL),0ULL);
-            u *= 0x1.0000000000000p-105L;
-        } else
-            SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL)-0x0690000000000000LL,0ULL);
+            SET_LDOUBLE_WORDS64(u, (hx & 0x7ff0000000000000LL), 0ULL);
+            u *= eps;
+        }
+        else
+            SET_LDOUBLE_WORDS64(
+                    u, (hx & 0x7ff0000000000000LL) - 0x0690000000000000LL,
+                    0ULL);
         return x - u;
-    } else {                /* p >= 0, x += ulp */
-        if((hx==0x7fefffffffffffffLL)&&(lx==0x7c8ffffffffffffeLL))
-            return x+x; /* overflow, return +inf */
-        if ((npy_uint64) hx >= 0xfff0000000000000ULL) {
-            SET_LDOUBLE_WORDS64(u,0xffefffffffffffffLL,0xfc8ffffffffffffeLL);
+    }
+    else { /* p >= 0, x += ulp */
+        if ((hx == 0x7fefffffffffffffLL) && (lx == 0x7c8ffffffffffffeLL))
+            return x + x; /* overflow, return +inf */
+        if ((npy_uint64)hx >= 0xfff0000000000000ULL) {
+            SET_LDOUBLE_WORDS64(u, 0xffefffffffffffffLL, 0xfc8ffffffffffffeLL);
             return u;
         }
-        if(ihx <= 0x0360000000000000LL) {  /* x <= LDBL_MIN */
-            u = math_opt_barrier (x);
+        if (ihx <= 0x0360000000000000LL) { /* x <= LDBL_MIN */
+            u = math_opt_barrier(x);
             x += LDBL_TRUE_MIN;
-            if (ihx < 0x0360000000000000LL
-                    || (hx > 0 && (npy_int64) lx < 0 && lx != 0x8000000000000001LL)
-                    || (hx < 0 && (npy_int64) lx >= 0)) {
+            if (ihx < 0x0360000000000000LL ||
+                (hx > 0 && (npy_int64)lx < 0 && lx != 0x8000000000000001LL) ||
+                (hx < 0 && (npy_int64)lx >= 0)) {
                 u = u * u;
-                math_force_eval (u);        /* raise underflow flag */
+                math_force_eval(u); /* raise underflow flag */
             }
-            if (x == 0.0L)  /* handle negative LDBL_TRUE_MIN case */
+            if (x == 0.0L) /* handle negative LDBL_TRUE_MIN case */
                 x = -0.0L;
             return x;
         }
         if (ihx < 0x06a0000000000000LL) { /* ulp will denormal */
-            SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL),0ULL);
-            u *= 0x1.0000000000000p-105L;
-        } else
-            SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL)-0x0690000000000000LL,0ULL);
+            SET_LDOUBLE_WORDS64(u, (hx & 0x7ff0000000000000LL), 0ULL);
+            u *= eps;
+        }
+        else
+            SET_LDOUBLE_WORDS64(
+                    u, (hx & 0x7ff0000000000000LL) - 0x0690000000000000LL,
+                    0ULL);
         return x + u;
     }
 }
 #else
-static npy_longdouble _nextl(npy_longdouble x, int p)
+static long double
+_next(long double x, int p)
 {
     volatile npy_longdouble t;
     union IEEEl2bitsrep ux;
@@ -266,53 +288,54 @@ static npy_longdouble _nextl(npy_longdouble x, int p)
 
     if ((GET_LDOUBLE_EXP(ux) == 0x7fff &&
          ((GET_LDOUBLE_MANH(ux) & ~LDBL_NBIT) | GET_LDOUBLE_MANL(ux)) != 0)) {
-        return ux.e;        /* x is nan */
+        return ux.e; /* x is nan */
     }
     if (ux.e == 0.0) {
-        SET_LDOUBLE_MANH(ux, 0);              /* return +-minsubnormal */
+        SET_LDOUBLE_MANH(ux, 0); /* return +-minsubnormal */
         SET_LDOUBLE_MANL(ux, 1);
         if (p >= 0) {
             SET_LDOUBLE_SIGN(ux, 0);
-        } else {
+        }
+        else {
             SET_LDOUBLE_SIGN(ux, 1);
         }
         t = ux.e * ux.e;
         if (t == ux.e) {
             return t;
-        } else {
-            return ux.e;           /* raise underflow flag */
+        }
+        else {
+            return ux.e; /* raise underflow flag */
         }
     }
-    if (p < 0) {      /* x -= ulp */
+    if (p < 0) { /* x -= ulp */
         if (GET_LDOUBLE_MANL(ux) == 0) {
             if ((GET_LDOUBLE_MANH(ux) & ~LDBL_NBIT) == 0) {
                 SET_LDOUBLE_EXP(ux, GET_LDOUBLE_EXP(ux) - 1);
             }
-            SET_LDOUBLE_MANH(ux,
-                             (GET_LDOUBLE_MANH(ux) - 1) |
-                             (GET_LDOUBLE_MANH(ux) & LDBL_NBIT));
+            SET_LDOUBLE_MANH(ux, (GET_LDOUBLE_MANH(ux) - 1) |
+                                         (GET_LDOUBLE_MANH(ux) & LDBL_NBIT));
         }
         SET_LDOUBLE_MANL(ux, GET_LDOUBLE_MANL(ux) - 1);
-    } else {                    /* x += ulp */
+    }
+    else { /* x += ulp */
         SET_LDOUBLE_MANL(ux, GET_LDOUBLE_MANL(ux) + 1);
         if (GET_LDOUBLE_MANL(ux) == 0) {
-            SET_LDOUBLE_MANH(ux,
-                             (GET_LDOUBLE_MANH(ux) + 1) |
-                             (GET_LDOUBLE_MANH(ux) & LDBL_NBIT));
+            SET_LDOUBLE_MANH(ux, (GET_LDOUBLE_MANH(ux) + 1) |
+                                         (GET_LDOUBLE_MANH(ux) & LDBL_NBIT));
             if ((GET_LDOUBLE_MANH(ux) & ~LDBL_NBIT) == 0) {
                 SET_LDOUBLE_EXP(ux, GET_LDOUBLE_EXP(ux) + 1);
             }
         }
     }
     if (GET_LDOUBLE_EXP(ux) == 0x7fff) {
-        return ux.e + ux.e;           /* overflow  */
+        return ux.e + ux.e; /* overflow  */
     }
-    if (GET_LDOUBLE_EXP(ux) == 0) {            /* underflow */
+    if (GET_LDOUBLE_EXP(ux) == 0) { /* underflow */
         if (LDBL_NBIT) {
             SET_LDOUBLE_MANH(ux, GET_LDOUBLE_MANH(ux) & ~LDBL_NBIT);
         }
         t = ux.e * ux.e;
-        if (t != ux.e) {           /* raise underflow flag */
+        if (t != ux.e) { /* raise underflow flag */
             return ux.e;
         }
     }
@@ -336,7 +359,8 @@ static npy_longdouble _nextl(npy_longdouble x, int p)
  */
 
 #ifndef HAVE_NEXTAFTER
-double npy_nextafter(double x, double y)
+double
+npy_nextafter(double x, double y)
 {
     volatile double t;
     npy_int32 hx, hy, ix, iy;
@@ -344,38 +368,42 @@ double npy_nextafter(double x, double y)
 
     EXTRACT_WORDS(hx, lx, x);
     EXTRACT_WORDS(hy, ly, y);
-    ix = hx & 0x7fffffff;       /* |x| */
-    iy = hy & 0x7fffffff;       /* |y| */
+    ix = hx & 0x7fffffff; /* |x| */
+    iy = hy & 0x7fffffff; /* |y| */
 
-    if (((ix >= 0x7ff00000) && ((ix - 0x7ff00000) | lx) != 0) ||        /* x is nan */
-        ((iy >= 0x7ff00000) && ((iy - 0x7ff00000) | ly) != 0))  /* y is nan */
+    if (((ix >= 0x7ff00000) && ((ix - 0x7ff00000) | lx) != 0) || /* x is nan */
+        ((iy >= 0x7ff00000) && ((iy - 0x7ff00000) | ly) != 0))   /* y is nan */
         return x + y;
     if (x == y)
-        return y;               /* x=y, return y */
-    if ((ix | lx) == 0) {       /* x == 0 */
-        INSERT_WORDS(x, hy & 0x80000000, 1);    /* return +-minsubnormal */
+        return y;                            /* x=y, return y */
+    if ((ix | lx) == 0) {                    /* x == 0 */
+        INSERT_WORDS(x, hy & 0x80000000, 1); /* return +-minsubnormal */
         t = x * x;
         if (t == x)
             return t;
         else
-            return x;           /* raise underflow flag */
+            return x; /* raise underflow flag */
     }
-    if (hx >= 0) {              /* x > 0 */
-        if (hx > hy || ((hx == hy) && (lx > ly))) {     /* x > y, x -= ulp */
+    if (hx >= 0) {                                  /* x > 0 */
+        if (hx > hy || ((hx == hy) && (lx > ly))) { /* x > y, x -= ulp */
             if (lx == 0)
                 hx -= 1;
             lx -= 1;
-        } else {                /* x < y, x += ulp */
+        }
+        else { /* x < y, x += ulp */
             lx += 1;
             if (lx == 0)
                 hx += 1;
         }
-    } else {                    /* x < 0 */
-        if (hy >= 0 || hx > hy || ((hx == hy) && (lx > ly))) {  /* x < y, x -= ulp */
+    }
+    else { /* x < 0 */
+        if (hy >= 0 || hx > hy ||
+            ((hx == hy) && (lx > ly))) { /* x < y, x -= ulp */
             if (lx == 0)
                 hx -= 1;
             lx -= 1;
-        } else {                /* x > y, x += ulp */
+        }
+        else { /* x > y, x += ulp */
             lx += 1;
             if (lx == 0)
                 hx += 1;
@@ -383,10 +411,10 @@ double npy_nextafter(double x, double y)
     }
     hy = hx & 0x7ff00000;
     if (hy >= 0x7ff00000)
-        return x + x;           /* overflow  */
-    if (hy < 0x00100000) {      /* underflow */
+        return x + x;      /* overflow  */
+    if (hy < 0x00100000) { /* underflow */
         t = x * x;
-        if (t != x) {           /* raise underflow flag */
+        if (t != x) { /* raise underflow flag */
             INSERT_WORDS(y, hx, lx);
             return y;
         }
@@ -397,48 +425,52 @@ double npy_nextafter(double x, double y)
 #endif
 
 #ifndef HAVE_NEXTAFTERF
-float npy_nextafterf(float x, float y)
+float
+npy_nextafterf(float x, float y)
 {
     volatile float t;
     npy_int32 hx, hy, ix, iy;
 
     GET_FLOAT_WORD(hx, x);
     GET_FLOAT_WORD(hy, y);
-    ix = hx & 0x7fffffff;       /* |x| */
-    iy = hy & 0x7fffffff;       /* |y| */
+    ix = hx & 0x7fffffff; /* |x| */
+    iy = hy & 0x7fffffff; /* |y| */
 
-    if ((ix > 0x7f800000) ||    /* x is nan */
-        (iy > 0x7f800000))      /* y is nan */
+    if ((ix > 0x7f800000) || /* x is nan */
+        (iy > 0x7f800000))   /* y is nan */
         return x + y;
     if (x == y)
-        return y;               /* x=y, return y */
-    if (ix == 0) {              /* x == 0 */
+        return y;                                 /* x=y, return y */
+    if (ix == 0) {                                /* x == 0 */
         SET_FLOAT_WORD(x, (hy & 0x80000000) | 1); /* return +-minsubnormal */
         t = x * x;
         if (t == x)
             return t;
         else
-            return x;           /* raise underflow flag */
+            return x; /* raise underflow flag */
     }
-    if (hx >= 0) {              /* x > 0 */
-        if (hx > hy) {          /* x > y, x -= ulp */
+    if (hx >= 0) {     /* x > 0 */
+        if (hx > hy) { /* x > y, x -= ulp */
             hx -= 1;
-        } else {                /* x < y, x += ulp */
+        }
+        else { /* x < y, x += ulp */
             hx += 1;
         }
-    } else {                    /* x < 0 */
-        if (hy >= 0 || hx > hy) {       /* x < y, x -= ulp */
+    }
+    else {                        /* x < 0 */
+        if (hy >= 0 || hx > hy) { /* x < y, x -= ulp */
             hx -= 1;
-        } else {                /* x > y, x += ulp */
+        }
+        else { /* x > y, x += ulp */
             hx += 1;
         }
     }
     hy = hx & 0x7f800000;
     if (hy >= 0x7f800000)
-        return x + x;           /* overflow  */
-    if (hy < 0x00800000) {      /* underflow */
+        return x + x;      /* overflow  */
+    if (hy < 0x00800000) { /* underflow */
         t = x * x;
-        if (t != x) {           /* raise underflow flag */
+        if (t != x) { /* raise underflow flag */
             SET_FLOAT_WORD(y, hx);
             return y;
         }
@@ -449,7 +481,8 @@ float npy_nextafterf(float x, float y)
 #endif
 
 #ifndef HAVE_NEXTAFTERL
-npy_longdouble npy_nextafterl(npy_longdouble x, npy_longdouble y)
+npy_longdouble
+npy_nextafterl(npy_longdouble x, npy_longdouble y)
 {
     volatile npy_longdouble t;
     union IEEEl2bitsrep ux;
@@ -462,52 +495,52 @@ npy_longdouble npy_nextafterl(npy_longdouble x, npy_longdouble y)
          ((GET_LDOUBLE_MANH(ux) & ~LDBL_NBIT) | GET_LDOUBLE_MANL(ux)) != 0) ||
         (GET_LDOUBLE_EXP(uy) == 0x7fff &&
          ((GET_LDOUBLE_MANH(uy) & ~LDBL_NBIT) | GET_LDOUBLE_MANL(uy)) != 0)) {
-        return ux.e + uy.e;        /* x or y is nan */
+        return ux.e + uy.e; /* x or y is nan */
     }
     if (ux.e == uy.e) {
-        return uy.e;               /* x=y, return y */
+        return uy.e; /* x=y, return y */
     }
     if (ux.e == 0.0) {
-        SET_LDOUBLE_MANH(ux, 0);              /* return +-minsubnormal */
+        SET_LDOUBLE_MANH(ux, 0); /* return +-minsubnormal */
         SET_LDOUBLE_MANL(ux, 1);
         SET_LDOUBLE_SIGN(ux, GET_LDOUBLE_SIGN(uy));
         t = ux.e * ux.e;
         if (t == ux.e) {
             return t;
-        } else {
-            return ux.e;           /* raise underflow flag */
+        }
+        else {
+            return ux.e; /* raise underflow flag */
         }
     }
-    if ((ux.e > 0.0) ^ (ux.e < uy.e)) {      /* x -= ulp */
+    if ((ux.e > 0.0) ^ (ux.e < uy.e)) { /* x -= ulp */
         if (GET_LDOUBLE_MANL(ux) == 0) {
             if ((GET_LDOUBLE_MANH(ux) & ~LDBL_NBIT) == 0) {
                 SET_LDOUBLE_EXP(ux, GET_LDOUBLE_EXP(ux) - 1);
             }
-            SET_LDOUBLE_MANH(ux,
-                             (GET_LDOUBLE_MANH(ux) - 1) |
-                             (GET_LDOUBLE_MANH(ux) & LDBL_NBIT));
+            SET_LDOUBLE_MANH(ux, (GET_LDOUBLE_MANH(ux) - 1) |
+                                         (GET_LDOUBLE_MANH(ux) & LDBL_NBIT));
         }
         SET_LDOUBLE_MANL(ux, GET_LDOUBLE_MANL(ux) - 1);
-    } else {                    /* x += ulp */
+    }
+    else { /* x += ulp */
         SET_LDOUBLE_MANL(ux, GET_LDOUBLE_MANL(ux) + 1);
         if (GET_LDOUBLE_MANL(ux) == 0) {
-            SET_LDOUBLE_MANH(ux,
-                             (GET_LDOUBLE_MANH(ux) + 1) |
-                             (GET_LDOUBLE_MANH(ux) & LDBL_NBIT));
+            SET_LDOUBLE_MANH(ux, (GET_LDOUBLE_MANH(ux) + 1) |
+                                         (GET_LDOUBLE_MANH(ux) & LDBL_NBIT));
             if ((GET_LDOUBLE_MANH(ux) & ~LDBL_NBIT) == 0) {
                 SET_LDOUBLE_EXP(ux, GET_LDOUBLE_EXP(ux) + 1);
             }
         }
     }
     if (GET_LDOUBLE_EXP(ux) == 0x7fff) {
-        return ux.e + ux.e;           /* overflow  */
+        return ux.e + ux.e; /* overflow  */
     }
-    if (GET_LDOUBLE_EXP(ux) == 0) {            /* underflow */
+    if (GET_LDOUBLE_EXP(ux) == 0) { /* underflow */
         if (LDBL_NBIT) {
             SET_LDOUBLE_MANH(ux, GET_LDOUBLE_MANH(ux) & ~LDBL_NBIT);
         }
         t = ux.e * ux.e;
-        if (t != ux.e) {           /* raise underflow flag */
+        if (t != ux.e) { /* raise underflow flag */
             return ux.e;
         }
     }
@@ -516,53 +549,100 @@ npy_longdouble npy_nextafterl(npy_longdouble x, npy_longdouble y)
 }
 #endif
 
-/**begin repeat
- * #suff = f,,l#
- * #SUFF = F,,L#
- * #type = npy_float, npy_double, npy_longdouble#
- */
-@type@ npy_spacing@suff@(@type@ x)
+namespace {
+template <typename T>
+struct numeric_limits;
+
+template <>
+struct numeric_limits<float> {
+    static const npy_float nan;
+};
+const npy_float numeric_limits<float>::nan = NPY_NANF;
+
+template <>
+struct numeric_limits<double> {
+    static const npy_double nan;
+};
+const npy_double numeric_limits<double>::nan = NPY_NAN;
+
+template <>
+struct numeric_limits<long double> {
+    static const npy_longdouble nan;
+};
+const npy_longdouble numeric_limits<long double>::nan = NPY_NANL;
+}  // namespace
+
+template <typename type>
+static type
+_npy_spacing(type x)
 {
     /* XXX: npy isnan/isinf may be optimized by bit twiddling */
     if (npy_isinf(x)) {
-        return NPY_NAN@SUFF@;
+        return numeric_limits<type>::nan;
     }
 
-    return _next@suff@(x, 1) - x;
+    return _next(x, 1) - x;
+}
+
+/*
+ * Instantiation of C interface
+ */
+extern "C" {
+npy_float
+npy_spacingf(npy_float x)
+{
+    return _npy_spacing(x);
+}
+npy_double
+npy_spacing(npy_double x)
+{
+    return _npy_spacing(x);
+}
+npy_longdouble
+npy_spacingl(npy_longdouble x)
+{
+    return _npy_spacing(x);
+}
 }
-/**end repeat**/
 
 /*
  * Decorate all the math functions which are available on the current platform
  */
 
 #ifdef HAVE_NEXTAFTERF
-float npy_nextafterf(float x, float y)
+extern "C" float
+npy_nextafterf(float x, float y)
 {
     return nextafterf(x, y);
 }
 #endif
 
 #ifdef HAVE_NEXTAFTER
-double npy_nextafter(double x, double y)
+extern "C" double
+npy_nextafter(double x, double y)
 {
     return nextafter(x, y);
 }
 #endif
 
 #ifdef HAVE_NEXTAFTERL
-npy_longdouble npy_nextafterl(npy_longdouble x, npy_longdouble y)
+extern "C" npy_longdouble
+npy_nextafterl(npy_longdouble x, npy_longdouble y)
 {
     return nextafterl(x, y);
 }
 #endif
 
-int npy_clear_floatstatus() {
-    char x=0;
+extern "C" int
+npy_clear_floatstatus()
+{
+    char x = 0;
     return npy_clear_floatstatus_barrier(&x);
 }
-int npy_get_floatstatus() {
-    char x=0;
+extern "C" int
+npy_get_floatstatus()
+{
+    char x = 0;
     return npy_get_floatstatus_barrier(&x);
 }
 
@@ -574,7 +654,6 @@ int npy_get_floatstatus() {
 #include <sys/param.h>
 #endif
 
-
 /*
  * Define floating point status functions. We must define
  * npy_get_floatstatus_barrier, npy_clear_floatstatus_barrier,
@@ -582,31 +661,32 @@ int npy_get_floatstatus() {
  * for all supported platforms.
  */
 
-
 /* Solaris --------------------------------------------------------*/
 /* --------ignoring SunOS ieee_flags approach, someone else can
 **         deal with that! */
-#if defined(sun) || defined(__BSD__) || defined(__OpenBSD__) || \
-    (defined(__FreeBSD__) && (__FreeBSD_version < 502114)) || \
-    defined(__NetBSD__)
+#if defined(sun) || defined(__BSD__) || defined(__OpenBSD__) ||   \
+        (defined(__FreeBSD__) && (__FreeBSD_version < 502114)) || \
+        defined(__NetBSD__)
 #include <ieeefp.h>
 
-int npy_get_floatstatus_barrier(char * param)
+extern "C" int
+npy_get_floatstatus_barrier(char *param)
 {
     int fpstatus = fpgetsticky();
     /*
      * By using a volatile, the compiler cannot reorder this call
      */
     if (param != NULL) {
-        volatile char NPY_UNUSED(c) = *(char*)param;
+        volatile char NPY_UNUSED(c) = *(char *)param;
     }
-    return ((FP_X_DZ  & fpstatus) ? NPY_FPE_DIVIDEBYZERO : 0) |
+    return ((FP_X_DZ & fpstatus) ? NPY_FPE_DIVIDEBYZERO : 0) |
            ((FP_X_OFL & fpstatus) ? NPY_FPE_OVERFLOW : 0) |
            ((FP_X_UFL & fpstatus) ? NPY_FPE_UNDERFLOW : 0) |
            ((FP_X_INV & fpstatus) ? NPY_FPE_INVALID : 0);
 }
 
-int npy_clear_floatstatus_barrier(char * param)
+extern "C" int
+npy_clear_floatstatus_barrier(char *param)
 {
     int fpstatus = npy_get_floatstatus_barrier(param);
     fpsetsticky(0);
@@ -614,38 +694,43 @@ int npy_clear_floatstatus_barrier(char * param)
     return fpstatus;
 }
 
-void npy_set_floatstatus_divbyzero(void)
+extern "C" void
+npy_set_floatstatus_divbyzero(void)
 {
     fpsetsticky(FP_X_DZ);
 }
 
-void npy_set_floatstatus_overflow(void)
+extern "C" void
+npy_set_floatstatus_overflow(void)
 {
     fpsetsticky(FP_X_OFL);
 }
 
-void npy_set_floatstatus_underflow(void)
+extern "C" void
+npy_set_floatstatus_underflow(void)
 {
     fpsetsticky(FP_X_UFL);
 }
 
-void npy_set_floatstatus_invalid(void)
+extern "C" void
+npy_set_floatstatus_invalid(void)
 {
     fpsetsticky(FP_X_INV);
 }
 
 #elif defined(_AIX) && !defined(__GNUC__)
-#include <float.h>
+#include <cfloat>
 #include <fpxcp.h>
 
-int npy_get_floatstatus_barrier(char *param)
+extern "C" int
+npy_get_floatstatus_barrier(char *param)
 {
     int fpstatus = fp_read_flag();
     /*
      * By using a volatile, the compiler cannot reorder this call
      */
     if (param != NULL) {
-        volatile char NPY_UNUSED(c) = *(char*)param;
+        volatile char NPY_UNUSED(c) = *(char *)param;
     }
     return ((FP_DIV_BY_ZERO & fpstatus) ? NPY_FPE_DIVIDEBYZERO : 0) |
            ((FP_OVERFLOW & fpstatus) ? NPY_FPE_OVERFLOW : 0) |
@@ -653,7 +738,8 @@ int npy_get_floatstatus_barrier(char *param)
            ((FP_INVALID & fpstatus) ? NPY_FPE_INVALID : 0);
 }
 
-int npy_clear_floatstatus_barrier(char * param)
+extern "C" int
+npy_clear_floatstatus_barrier(char *param)
 {
     int fpstatus = npy_get_floatstatus_barrier(param);
     fp_swap_flag(0);
@@ -661,28 +747,32 @@ int npy_clear_floatstatus_barrier(char * param)
     return fpstatus;
 }
 
-void npy_set_floatstatus_divbyzero(void)
+extern "C" void
+npy_set_floatstatus_divbyzero(void)
 {
     fp_raise_xcp(FP_DIV_BY_ZERO);
 }
 
-void npy_set_floatstatus_overflow(void)
+extern "C" void
+npy_set_floatstatus_overflow(void)
 {
     fp_raise_xcp(FP_OVERFLOW);
 }
 
-void npy_set_floatstatus_underflow(void)
+extern "C" void
+npy_set_floatstatus_underflow(void)
 {
     fp_raise_xcp(FP_UNDERFLOW);
 }
 
-void npy_set_floatstatus_invalid(void)
+extern "C" void
+npy_set_floatstatus_invalid(void)
 {
     fp_raise_xcp(FP_INVALID);
 }
 
 #elif defined(_MSC_VER) || (defined(__osf__) && defined(__alpha)) || \
-      defined (__UCLIBC__) || (defined(__arc__) && defined(__GLIBC__))
+        defined(__UCLIBC__) || (defined(__arc__) && defined(__GLIBC__))
 
 /*
  * By using a volatile floating point value,
@@ -694,25 +784,29 @@ void npy_set_floatstatus_invalid(void)
  * a race condition.
  */
 static volatile double _npy_floatstatus_x,
-    _npy_floatstatus_zero = 0.0, _npy_floatstatus_big = 1e300,
-    _npy_floatstatus_small = 1e-300, _npy_floatstatus_inf;
+        _npy_floatstatus_zero = 0.0, _npy_floatstatus_big = 1e300,
+        _npy_floatstatus_small = 1e-300, _npy_floatstatus_inf;
 
-void npy_set_floatstatus_divbyzero(void)
+extern "C" void
+npy_set_floatstatus_divbyzero(void)
 {
     _npy_floatstatus_x = 1.0 / _npy_floatstatus_zero;
 }
 
-void npy_set_floatstatus_overflow(void)
+extern "C" void
+npy_set_floatstatus_overflow(void)
 {
     _npy_floatstatus_x = _npy_floatstatus_big * 1e300;
 }
 
-void npy_set_floatstatus_underflow(void)
+extern "C" void
+npy_set_floatstatus_underflow(void)
 {
     _npy_floatstatus_x = _npy_floatstatus_small * 1e-300;
 }
 
-void npy_set_floatstatus_invalid(void)
+extern "C" void
+npy_set_floatstatus_invalid(void)
 {
     _npy_floatstatus_inf = NPY_INFINITY;
     _npy_floatstatus_x = _npy_floatstatus_inf - NPY_INFINITY;
@@ -721,9 +815,10 @@ void npy_set_floatstatus_invalid(void)
 /* MS Windows -----------------------------------------------------*/
 #if defined(_MSC_VER)
 
-#include <float.h>
+#include <cfloat>
 
-int npy_get_floatstatus_barrier(char *param)
+extern "C" int
+npy_get_floatstatus_barrier(char *param)
 {
     /*
      * By using a volatile, the compiler cannot reorder this call
@@ -732,12 +827,12 @@ int npy_get_floatstatus_barrier(char *param)
     int fpstatus = _statusfp();
 #else
     /* windows enables sse on 32 bit, so check both flags */
-    int fpstatus, fpstatus2;
+    unsigned int fpstatus, fpstatus2;
     _statusfp2(&fpstatus, &fpstatus2);
     fpstatus |= fpstatus2;
 #endif
     if (param != NULL) {
-        volatile char NPY_UNUSED(c) = *(char*)param;
+        volatile char NPY_UNUSED(c) = *(char *)param;
     }
     return ((SW_ZERODIVIDE & fpstatus) ? NPY_FPE_DIVIDEBYZERO : 0) |
            ((SW_OVERFLOW & fpstatus) ? NPY_FPE_OVERFLOW : 0) |
@@ -745,7 +840,8 @@ int npy_get_floatstatus_barrier(char *param)
            ((SW_INVALID & fpstatus) ? NPY_FPE_INVALID : 0);
 }
 
-int npy_clear_floatstatus_barrier(char *param)
+extern "C" int
+npy_clear_floatstatus_barrier(char *param)
 {
     int fpstatus = npy_get_floatstatus_barrier(param);
     _clearfp();
@@ -758,22 +854,24 @@ int npy_clear_floatstatus_barrier(char *param)
 
 #include <machine/fpu.h>
 
-int npy_get_floatstatus_barrier(char *param)
+extern "C" int
+npy_get_floatstatus_barrier(char *param)
 {
     unsigned long fpstatus = ieee_get_fp_control();
     /*
      * By using a volatile, the compiler cannot reorder this call
      */
     if (param != NULL) {
-        volatile char NPY_UNUSED(c) = *(char*)param;
+        volatile char NPY_UNUSED(c) = *(char *)param;
     }
-    return  ((IEEE_STATUS_DZE & fpstatus) ? NPY_FPE_DIVIDEBYZERO : 0) |
-            ((IEEE_STATUS_OVF & fpstatus) ? NPY_FPE_OVERFLOW : 0) |
-            ((IEEE_STATUS_UNF & fpstatus) ? NPY_FPE_UNDERFLOW : 0) |
-            ((IEEE_STATUS_INV & fpstatus) ? NPY_FPE_INVALID : 0);
+    return ((IEEE_STATUS_DZE & fpstatus) ? NPY_FPE_DIVIDEBYZERO : 0) |
+           ((IEEE_STATUS_OVF & fpstatus) ? NPY_FPE_OVERFLOW : 0) |
+           ((IEEE_STATUS_UNF & fpstatus) ? NPY_FPE_UNDERFLOW : 0) |
+           ((IEEE_STATUS_INV & fpstatus) ? NPY_FPE_INVALID : 0);
 }
 
-int npy_clear_floatstatus_barrier(char *param)
+extern "C" int
+npy_clear_floatstatus_barrier(char *param)
 {
     int fpstatus = npy_get_floatstatus_barrier(param);
     /* clear status bits as well as disable exception mode if on */
@@ -787,54 +885,58 @@ int npy_clear_floatstatus_barrier(char *param)
 
 #else
 /* General GCC code, should work on most platforms */
-#  include <fenv.h>
+#include <fenv.h>
 
-int npy_get_floatstatus_barrier(char* param)
+extern "C" int
+npy_get_floatstatus_barrier(char *param)
 {
-    int fpstatus = fetestexcept(FE_DIVBYZERO | FE_OVERFLOW |
-                                FE_UNDERFLOW | FE_INVALID);
+    int fpstatus = fetestexcept(FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW |
+                                FE_INVALID);
     /*
      * By using a volatile, the compiler cannot reorder this call
      */
     if (param != NULL) {
-        volatile char NPY_UNUSED(c) = *(char*)param;
+        volatile char NPY_UNUSED(c) = *(char *)param;
     }
 
-    return ((FE_DIVBYZERO  & fpstatus) ? NPY_FPE_DIVIDEBYZERO : 0) |
-           ((FE_OVERFLOW   & fpstatus) ? NPY_FPE_OVERFLOW : 0) |
-           ((FE_UNDERFLOW  & fpstatus) ? NPY_FPE_UNDERFLOW : 0) |
-           ((FE_INVALID    & fpstatus) ? NPY_FPE_INVALID : 0);
+    return ((FE_DIVBYZERO & fpstatus) ? NPY_FPE_DIVIDEBYZERO : 0) |
+           ((FE_OVERFLOW & fpstatus) ? NPY_FPE_OVERFLOW : 0) |
+           ((FE_UNDERFLOW & fpstatus) ? NPY_FPE_UNDERFLOW : 0) |
+           ((FE_INVALID & fpstatus) ? NPY_FPE_INVALID : 0);
 }
 
-int npy_clear_floatstatus_barrier(char * param)
+extern "C" int
+npy_clear_floatstatus_barrier(char *param)
 {
     /* testing float status is 50-100 times faster than clearing on x86 */
     int fpstatus = npy_get_floatstatus_barrier(param);
     if (fpstatus != 0) {
-        feclearexcept(FE_DIVBYZERO | FE_OVERFLOW |
-                      FE_UNDERFLOW | FE_INVALID);
+        feclearexcept(FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW | FE_INVALID);
     }
 
     return fpstatus;
 }
 
-
-void npy_set_floatstatus_divbyzero(void)
+extern "C" void
+npy_set_floatstatus_divbyzero(void)
 {
     feraiseexcept(FE_DIVBYZERO);
 }
 
-void npy_set_floatstatus_overflow(void)
+extern "C" void
+npy_set_floatstatus_overflow(void)
 {
     feraiseexcept(FE_OVERFLOW);
 }
 
-void npy_set_floatstatus_underflow(void)
+extern "C" void
+npy_set_floatstatus_underflow(void)
 {
     feraiseexcept(FE_UNDERFLOW);
 }
 
-void npy_set_floatstatus_invalid(void)
+extern "C" void
+npy_set_floatstatus_invalid(void)
 {
     feraiseexcept(FE_INVALID);
 }