BUG: do not use bitfields for bit-twidling of long doubles.

Some compilers (Sun CC v12 for example) do not support bitfields with
more than 32 bits. We redesign the internal macros and union usages to
use bitmask instead where needed.

Not thoroughly tested yet.

2 files changed, 233 insertions, 99 deletions

diff --git a/numpy/core/src/npymath/ieee754.c.src b/numpy/core/src/npymath/ieee754.c.src
index 52a7923da..6a27ee58c 100644
--- a/numpy/core/src/npymath/ieee754.c.src
+++ b/numpy/core/src/npymath/ieee754.c.src
@@ -159,59 +159,66 @@ float npy_nextafterf(float x, float y)
 npy_longdouble npy_nextafterl(npy_longdouble x, npy_longdouble y)
 {
     volatile npy_longdouble t;
-    union IEEEl2bits ux, uy;
+    ldouble_man_t xmanl, xmanh, ymanl, ymanh;
+    ldouble_exp_t xexp, yexp;
+    ldouble_sign_t xsign, ysign;
 
-    ux.e = x;
-    uy.e = y;
+	GET_LDOUBLE_SIGN(x, xsign);
+	GET_LDOUBLE_EXP(x, xexp);
+	GET_LDOUBLE_MANL(x, xmanl);
+	GET_LDOUBLE_MANH(x, xmanh);
 
-    if ((ux.bits.exp == 0x7fff &&
-         ((ux.bits.manh & ~LDBL_NBIT) | ux.bits.manl) != 0) ||
-        (uy.bits.exp == 0x7fff &&
-         ((uy.bits.manh & ~LDBL_NBIT) | uy.bits.manl) != 0)) {
+	GET_LDOUBLE_SIGN(y, ysign);
+	GET_LDOUBLE_EXP(y, yexp);
+	GET_LDOUBLE_MANL(y, ymanl);
+	GET_LDOUBLE_MANH(y, ymanh);
+
+    if ((xexp == 0x7fff && ((xmanh & ~LDBL_NBIT) | xmanl) != 0) ||
+        (yexp == 0x7fff && ((ymanh & ~LDBL_NBIT) | ymanl) != 0)) {
         return x + y;           /* x or y is nan */
     }
     if (x == y) {
         return y;               /* x=y, return y */
     }
     if (x == 0.0) {
-        ux.bits.manh = 0;       /* return +-minsubnormal */
-        ux.bits.manl = 1;
-        ux.bits.sign = uy.bits.sign;
-        t = ux.e * ux.e;
-        if (t == ux.e) {
+        xmanh = 0;              /* return +-minsubnormal */
+        xmanl = 1;
+        xsign = ysign;
+        t = x * x;
+        if (t == x) {
             return t;
         } else {
-            return ux.e;        /* raise underflow flag */
+            return x;           /* raise underflow flag */
         }
     }
     if (x > 0.0 ^ x < y) {      /* x -= ulp */
-        if (ux.bits.manl == 0) {
-            if ((ux.bits.manh & ~LDBL_NBIT) == 0) {
-                ux.bits.exp -= 1;
+        if (xmanl == 0) {
+            if ((xmanh & ~LDBL_NBIT) == 0) {
+                xexp -= 1;
             }
-            ux.bits.manh = (ux.bits.manh - 1) | (ux.bits.manh & LDBL_NBIT);
+            xmanh = (xmanh - 1) | (xmanh & LDBL_NBIT);
         }
-        ux.bits.manl -= 1;
+        xmanl -= 1;
     } else {                    /* x += ulp */
-        ux.bits.manl += 1;
-        if (ux.bits.manl == 0) {
-            ux.bits.manh = (ux.bits.manh + 1) | (ux.bits.manh & LDBL_NBIT);
-            if ((ux.bits.manh & ~LDBL_NBIT) == 0) {
-                ux.bits.exp += 1;
+        xmanl += 1;
+        if (xmanl == 0) {
+            xmanh = (xmanh + 1) | (xmanh & LDBL_NBIT);
+            if ((xmanh & ~LDBL_NBIT) == 0) {
+                xexp += 1;
             }
         }
     }
-    if (ux.bits.exp == 0x7fff) {
+    if (xexp == 0x7fff) {
         return x + x;           /* overflow  */
     }
-    if (ux.bits.exp == 0) {     /* underflow */
-        mask_nbit_l(ux);
-        t = ux.e * ux.e;
-        if (t != ux.e) {        /* raise underflow flag */
-            return ux.e;
+    if (xexp == 0) {            /* underflow */
+        mask_nbit_l(xmanh);
+        t = x * x;
+        if (t != x) {           /* raise underflow flag */
+            return x;
         }
     }
-    return ux.e;
+    return x;
 }
 #endif
 
diff --git a/numpy/core/src/npymath/npy_math_private.h b/numpy/core/src/npymath/npy_math_private.h
index 5e309a30b..9880b5ed8 100644
--- a/numpy/core/src/npymath/npy_math_private.h
+++ b/numpy/core/src/npymath/npy_math_private.h
@@ -175,88 +175,215 @@ do {								\
  * Long double support
  */
 #if defined(HAVE_LDOUBLE_INTEL_EXTENDED_12_BYTES_LE)
-    /* 80 bits Intel format aligned on 12 bytes */
-    union IEEEl2bits {
+    /*
+     * Intel extended 80 bits precision. Bit representation is
+     *          |  junk  |     s  |eeeeeeeeeeeeeee|mmmmmmmm................mmmmmmm|
+     *          | 16 bits|  1 bit |    15 bits    |            64 bits            |
+     *          |             a[2]                |     a[1]     |    a[0]        |
+     *
+     * 16 stronger bits of a[2] are junk
+     */
+    union IEEEl2bitsrep {
         npy_longdouble  e;
-        struct {
-            npy_uint32  manl    :32;
-            npy_uint32  manh    :32;
-            npy_uint32  exp     :15;
-            npy_uint32  sign    :1;
-            npy_uint32  junk    :16;
-        } bits;
+        npy_uint32      a[3];
     };
-    #define LDBL_NBIT   0x80000000
-    #define mask_nbit_l(u)  ((u).bits.manh &= ~LDBL_NBIT)
+
+    #define LDBL_MANL_INDEX     0
+    #define LDBL_MANL_MASK      0xFFFFFFFF
+    #define LDBL_MANL_SHIFT(x)  (x)
+
+    #define LDBL_MANH_INDEX     1
+    #define LDBL_MANH_MASK      0xFFFFFFFF
+    #define LDBL_MANH_SHIFT(x)  (x)
+
+    #define LDBL_EXP_INDEX      2
+    #define LDBL_EXP_MASK       0x7FFF
+    #define LDBL_EXP_SHIFT(x)   (x)
+
+    #define LDBL_SIGN_INDEX     2
+    #define LDBL_SIGN_MASK      0x8000
+    #define LDBL_SIGN_SHIFT(x)  ((x) >>  15)
+
+    #define LDBL_NBIT           0x80000000
+    #define mask_nbit_l(manh)   ((manh) &= ~LDBL_NBIT)
+
+    typedef npy_uint32 ldouble_man_t;
+    typedef npy_uint32 ldouble_exp_t;
+    typedef npy_uint32 ldouble_sign_t;
 #elif defined(HAVE_LDOUBLE_INTEL_EXTENDED_16_BYTES_LE)
-    /* 80 bits Intel format aligned on 16 bytes */
-    union IEEEl2bits {
-        npy_longdouble  e;
-        struct {
-            npy_uint32  manl    :32;
-            npy_uint32  manh    :32;
-            npy_uint32  exp     :15;
-            npy_uint32  sign    :1;
-            npy_uint32  junkl   :16;
-            npy_uint32  junkh   :32;
-        } bits;
+    /* 
+     * Intel extended 80 bits precision, 16 bytes alignment.. Bit representation is
+     *		|  junk	 |     s  |eeeeeeeeeeeeeee|mmmmmmmm................mmmmmmm|
+     *		| 16 bits|  1 bit |    15 bits    |            64 bits	          |
+     *		|             a[2]       	  |  	a[1]	 |    a[0]	  |
+     *
+     * a[3] and 16 stronger bits of a[2] are junk
+     */
+    union IEEEl2bitsrep {
+    	npy_longdouble e;
+    	npy_uint32     a[4];
     };
-    #define LDBL_NBIT   0x80000000
-    #define mask_nbit_l(u)  ((u).bits.manh &= ~LDBL_NBIT)
-#elif defined(HAVE_LDOUBLE_IEEE_DOUBLE_16_BYTES_BE)
+
+    #define LDBL_MANL_INDEX     0
+    #define LDBL_MANL_MASK 	0xFFFFFFFF
+    #define LDBL_MANL_SHIFT(x)  ((x))
+    
+    #define LDBL_MANH_INDEX     1
+    #define LDBL_MANH_MASK 	0xFFFFFFFF
+    #define LDBL_MANH_SHIFT(x)  ((x))
+    
+    #define LDBL_EXP_INDEX 	2
+    #define LDBL_EXP_MASK  	0x7FFF
+    #define LDBL_EXP_SHIFT(x)   (x)
+    
+    #define LDBL_SIGN_INDEX     2
+    #define LDBL_SIGN_MASK	0x8000
+    #define LDBL_SIGN_SHIFT(x)  ((x) >>  15)
+
+    #define LDBL_NBIT           0x80000000
+    #define mask_nbit_l(manh)   ((manh) &= ~LDBL_NBIT)
+
+    typedef npy_uint32 ldouble_man_t;
+    typedef npy_uint32 ldouble_exp_t;
+    typedef npy_uint32 ldouble_sign_t;
+#elif defined(HAVE_LDOUBLE_IEEE_DOUBLE_16_BYTES_BE) || \
+      defined(HAVE_LDOUBLE_IEEE_DOUBLE_BE)
     /* 64 bits IEEE double precision aligned on 16 bytes: used by ppc arch on
      * Mac OS X */
-    union IEEEl2bits {
-        npy_longdouble  e;
-        struct {
-            npy_uint32  sign    :1;
-            npy_uint32  exp     :11;
-            npy_uint32  manh    :20;
-            npy_uint32  manl    :32;
-            npy_uint32  junkh   :32;
-            npy_uint32  junkl   :32;
-        } bits;
-    };
-    #define LDBL_NBIT   0
-    #define mask_nbit_l(u)  ((void)0)
-#elif defined(HAVE_LDOUBLE_IEEE_QUAD_BE)
-    /* Quad precision IEEE format */
-    union IEEEl2bits {
-        npy_longdouble  e;
-        struct {
-            npy_uint32  sign    :1;
-            npy_uint32  exp     :15;
-            npy_uint64  manh    :48;
-            npy_uint64  manl    :64;
-        } bits;
+
+    /* 
+     * IEEE double precision. Bit representation is
+     *		|  s  |eeeeeeeeeee|mmmmmmmm................mmmmmmm|
+     *		|1 bit|  11 bits  |            52 bits		  |
+     *		|          a[0]         |  	  a[1]		  |
+     */
+    union IEEEl2bitsrep {
+    	long double e;
+    	npy_uint32 a[2];
     };
-    #define mask_nbit_l(u)  ((void)0)
+    
+    #define LDBL_MANL_INDEX     1
+    #define LDBL_MANL_MASK 	0xFFFFFFFF
+    #define LDBL_MANL_SHIFT(x)  ((x))
+    
+    #define LDBL_MANH_INDEX     0
+    #define LDBL_MANH_MASK 	0x000FFFFF
+    #define LDBL_MANH_SHIFT(x)  ((x))
+    
+    #define LDBL_EXP_INDEX 	0
+    #define LDBL_EXP_MASK  	0x7FF00000
+    #define LDBL_EXP_SHIFT(x)   ((x) >>  20)
+    
+    #define LDBL_SIGN_INDEX     0
+    #define LDBL_SIGN_MASK	0x80000000
+    #define LDBL_SIGN_SHIFT(x)  ((x) >>  31)
+
+    #define LDBL_NBIT           0
+    #define mask_nbit_l(u)      ((void)0)
+
+    typedef npy_uint32 ldouble_man_t;
+    typedef npy_uint32 ldouble_exp_t;
+    typedef npy_uint32 ldouble_sign_t;
 #elif defined(HAVE_LDOUBLE_IEEE_DOUBLE_LE)
-    union IEEEl2bits {
-        npy_longdouble  e;
-        struct {
-            npy_unit32  manl    :32;
-            npy_unit32  manh    :20;
-            npy_unit32  exp     :11;
-            npy_unit32  sign    :1;
-        } bits;
+    /* 64 bits IEEE double precision, Little Endian. */
+
+    /* 
+     * IEEE double precision. Bit representation is
+     *		|  s  |eeeeeeeeeee|mmmmmmmm................mmmmmmm|
+     *		|1 bit|  11 bits  |            52 bits		  |
+     *		|          a[1]         |  	  a[0]		  |
+     */
+    union IEEEl2bitsrep {
+    	npy_longdouble  e;
+    	npy_uint32      a[2];
     };
-    #define LDBL_NBIT   0
-    #define mask_nbit_l(u)  ((void)0)
-#elif defined(HAVE_LDOUBLE_IEEE_DOUBLE_BE)
-    union IEEEl2bits {
-        npy_longdouble  e;
-        struct {
-            npy_unit32  sign    :1;
-            npy_unit32  exp     :11;
-            npy_unit32  manh    :20;
-            npy_unit32  manl    :32;
-        } bits;
+    
+    #define LDBL_MANL_INDEX     0
+    #define LDBL_MANL_MASK 	0xFFFFFFFF
+    #define LDBL_MANL_SHIFT(x)  ((x))
+    
+    #define LDBL_MANH_INDEX     1
+    #define LDBL_MANH_MASK 	0x000FFFFF
+    #define LDBL_MANH_SHIFT(x)  ((x))
+    
+    #define LDBL_EXP_INDEX 	1
+    #define LDBL_EXP_MASK  	0x7FF00000
+    #define LDBL_EXP_SHIFT(x)   ((x) >>  20)
+    
+    #define LDBL_SIGN_INDEX     1
+    #define LDBL_SIGN_MASK	0x80000000
+    #define LDBL_SIGN_SHIFT(x)  ((x) >>  31)
+
+    #define LDBL_NBIT           0x80000000
+    #define mask_nbit_l(manh)   ((manh) &= ~LDBL_NBIT)
+
+    typedef npy_uint32 ldouble_man_t;
+    typedef npy_uint32 ldouble_exp_t;
+    typedef npy_uint32 ldouble_sign_t;
+#elif defined(HAVE_LDOUBLE_IEEE_QUAD_BE)
+    /* 
+     * IEEE quad precision, Big Endian. Bit representation is
+     *		|  s  |eeeeeeeeeee|mmmmmmmm................mmmmmmm|
+     *		|1 bit|  15 bits  |            112 bits		  |
+     *		|          a[0]         |           a[1]          |
+     */
+    union IEEEl2bitsrep {
+    	npy_longdouble  e;
+    	npy_uint64      a[2];
     };
-    #define LDBL_NBIT   0
+    
+    #define LDBL_MANL_INDEX     1
+    #define LDBL_MANL_MASK 	0xFFFFFFFFFFFFFFFF
+    #define LDBL_MANL_SHIFT(x)  ((x))
+    
+    #define LDBL_MANH_INDEX     0
+    #define LDBL_MANH_MASK 	0x0000FFFFFFFFFFFF
+    #define LDBL_MANH_SHIFT(x)  ((x))
+    
+    #define LDBL_EXP_INDEX 	0
+    #define LDBL_EXP_MASK  	0x7FFF000000000000
+    #define LDBL_EXP_SHIFT(x)   ((x) >>  48)
+    
+    #define LDBL_SIGN_INDEX     0
+    #define LDBL_SIGN_MASK	0x8000000000000000
+    #define LDBL_SIGN_SHIFT(x)  ((x) >>  63)
+
     #define mask_nbit_l(u)  ((void)0)
+
+    typedef npy_uint64 ldouble_man_t;
+    typedef npy_uint64 ldouble_exp_t;
+    typedef npy_uint32 ldouble_sign_t;
 #endif
 
+/* Put the sign bit of x into i - i should be defined as ldouble_sign_t */
+#define GET_LDOUBLE_SIGN(x, i)	\
+do {				\
+	const union IEEEl2bitsrep u = {x};	\
+	(i) = LDBL_SIGN_SHIFT(u.a[LDBL_SIGN_INDEX] & LDBL_SIGN_MASK);  \
+} while (0)
+
+/* Put the exp bits of x into i - i should be defined as ldouble_exp_t */
+#define GET_LDOUBLE_EXP(x, i)	\
+do {				\
+	const union IEEEl2bitsrep u = {x};	\
+	(i) = LDBL_EXP_SHIFT(u.a[LDBL_EXP_INDEX] & LDBL_EXP_MASK);  \
+} while (0)
+
+/* Put the manl bits of x into i - i should be defined as ldouble_man_t */
+#define GET_LDOUBLE_MANL(x, i)	\
+do {				\
+	const union IEEEl2bitsrep u = {x};	\
+	(i) = LDBL_MANL_SHIFT(u.a[LDBL_MANL_INDEX] & LDBL_MANL_MASK);  \
+} while (0)
+
+/* Put the manh bits (hist bits of mantisss) of x into i - i should be defined
+ * as ldouble_man_t */
+#define GET_LDOUBLE_MANH(x, i)	\
+do {				\
+	const union IEEEl2bitsrep u = {x};	\
+	(i) = LDBL_MANH_SHIFT(u.a[LDBL_MANH_INDEX] & LDBL_MANH_MASK);  \
+} while (0)
+
 /*
  * Those unions are used to convert a pointer of npy_cdouble to native C99
  * complex or our own complex type indenpendently on whether C99 complex