From 85ed88218164ceecb360bd9b8539e55076df9d19 Mon Sep 17 00:00:00 2001

From: =?UTF-8?q?Nikola=20Forr=C3=B3?= <nforro@redhat.com>

Date: Thu, 30 May 2019 15:52:23 +0200

Subject: [PATCH] Fix broken float128 on all arches except x86_64

Backport of the following upstream commits:

3d21112d6 MAINT: Add bitmask helper functions

4d7c529d7 MAINT: Reorganize Dragon4 code to clarify float formats

727756e63 MAINT: replace static BigInts in dragon4 by dummy manager

888d37086 MAINT: Add comments to long_double detection code

cdc6b6847 BUG: Implement float128 dragon4 for IBM double-double (ppc64)

718d9e866 ENH: Add support for the 64-bit RISC-V architecture

b35dfc2d5 ENH: Add AARCH32 support.

3973b2508 Adjust endianness header file to accommodate AARCHXX changes.

f07359b22 BLD: Modify cpu detection and function to get build working on aarch64  (#11568)

dd949d5c6 BUG: Fix printing of longdouble on ppc64le.

641058fb7 MAINT: remove darwin hardcoded LDOUBLE detection

ec843c700 BUG: Fix undefined functions on big-endian systems.

00bc5bc88 BUG: Fix test sensitive to platform byte order.

---

 numpy/core/include/numpy/npy_cpu.h          |   30 +-

 numpy/core/include/numpy/npy_endian.h       |   44 +-

 numpy/core/setup.py                         |   13 +-

 numpy/core/setup_common.py                  |   18 +-

 numpy/core/src/multiarray/dragon4.c         | 1613 ++++++++++++-------

 numpy/core/src/multiarray/dragon4.h         |   66 +-

 numpy/core/src/multiarray/scalartypes.c.src |    5 +-

 numpy/core/src/npymath/ieee754.c.src        |    3 +-

 numpy/core/src/npymath/npy_math_private.h   |    9 +-

 numpy/core/src/private/npy_config.h         |    3 +-

 numpy/core/src/private/npy_fpmath.h         |   41 +-

 numpy/core/tests/test_scalarprint.py        |   63 +-

 numpy/core/tests/test_umath.py              |   12 +-

 numpy/ma/tests/test_core.py                 |   11 +-

 14 files changed, 1277 insertions(+), 654 deletions(-)

diff --git a/numpy/core/include/numpy/npy_cpu.h b/numpy/core/include/numpy/npy_cpu.h

index 84653ea..1109426 100644

--- a/numpy/core/include/numpy/npy_cpu.h

+++ b/numpy/core/include/numpy/npy_cpu.h

@@ -17,6 +17,7 @@

  *              NPY_CPU_SH_BE

  *              NPY_CPU_ARCEL

  *              NPY_CPU_ARCEB

+ *              NPY_CPU_RISCV64

  */

 #ifndef _NPY_CPUARCH_H_

 #define _NPY_CPUARCH_H_

@@ -60,10 +61,27 @@

     #define NPY_CPU_HPPA

 #elif defined(__alpha__)

     #define NPY_CPU_ALPHA

-#elif defined(__arm__) && defined(__ARMEL__)

-    #define NPY_CPU_ARMEL

-#elif defined(__arm__) && defined(__ARMEB__)

-    #define NPY_CPU_ARMEB

+#elif defined(__arm__) || defined(__aarch64__)

+    #if defined(__ARMEB__) || defined(__AARCH64EB__)

+        #if defined(__ARM_32BIT_STATE)

+            #define NPY_CPU_ARMEB_AARCH32

+        #elif defined(__ARM_64BIT_STATE)

+            #define NPY_CPU_ARMEB_AARCH64

+        #else

+            #define NPY_CPU_ARMEB

+        #endif

+    #elif defined(__ARMEL__) || defined(__AARCH64EL__)

+        #if defined(__ARM_32BIT_STATE)

+            #define NPY_CPU_ARMEL_AARCH32

+        #elif defined(__ARM_64BIT_STATE)

+            #define NPY_CPU_ARMEL_AARCH64

+        #else

+            #define NPY_CPU_ARMEL

+        #endif

+    #else

+        # error Unknown ARM CPU, please report this to numpy maintainers with \

+	information about your platform (OS, CPU and compiler)

+    #endif

 #elif defined(__sh__) && defined(__LITTLE_ENDIAN__)

     #define NPY_CPU_SH_LE

 #elif defined(__sh__) && defined(__BIG_ENDIAN__)

@@ -74,14 +92,14 @@

     #define NPY_CPU_MIPSEB

 #elif defined(__or1k__)

     #define NPY_CPU_OR1K

-#elif defined(__aarch64__)

-    #define NPY_CPU_AARCH64

 #elif defined(__mc68000__)

     #define NPY_CPU_M68K

 #elif defined(__arc__) && defined(__LITTLE_ENDIAN__)

     #define NPY_CPU_ARCEL

 #elif defined(__arc__) && defined(__BIG_ENDIAN__)

     #define NPY_CPU_ARCEB

+#elif defined(__riscv) && defined(__riscv_xlen) && __riscv_xlen == 64

+    #define NPY_CPU_RISCV64

 #else

     #error Unknown CPU, please report this to numpy maintainers with \

     information about your platform (OS, CPU and compiler)

diff --git a/numpy/core/include/numpy/npy_endian.h b/numpy/core/include/numpy/npy_endian.h

index 1a42121..44cdffd 100644

--- a/numpy/core/include/numpy/npy_endian.h

+++ b/numpy/core/include/numpy/npy_endian.h

@@ -37,27 +37,31 @@

     #define NPY_LITTLE_ENDIAN 1234

     #define NPY_BIG_ENDIAN 4321

-    #if defined(NPY_CPU_X86)            \

-            || defined(NPY_CPU_AMD64)   \

-            || defined(NPY_CPU_IA64)    \

-            || defined(NPY_CPU_ALPHA)   \

-            || defined(NPY_CPU_ARMEL)   \

-            || defined(NPY_CPU_AARCH64) \

-            || defined(NPY_CPU_SH_LE)   \

-            || defined(NPY_CPU_MIPSEL)  \

-            || defined(NPY_CPU_PPC64LE) \

-            || defined(NPY_CPU_ARCEL)

+    #if defined(NPY_CPU_X86)                  \

+            || defined(NPY_CPU_AMD64)         \

+            || defined(NPY_CPU_IA64)          \

+            || defined(NPY_CPU_ALPHA)         \

+            || defined(NPY_CPU_ARMEL)         \

+            || defined(NPY_CPU_ARMEL_AARCH32) \

+            || defined(NPY_CPU_ARMEL_AARCH64) \

+            || defined(NPY_CPU_SH_LE)         \

+            || defined(NPY_CPU_MIPSEL)        \

+            || defined(NPY_CPU_PPC64LE)       \

+            || defined(NPY_CPU_ARCEL)         \

+            || defined(NPY_CPU_RISCV64)

         #define NPY_BYTE_ORDER NPY_LITTLE_ENDIAN

-    #elif defined(NPY_CPU_PPC)          \

-            || defined(NPY_CPU_SPARC)   \

-            || defined(NPY_CPU_S390)    \

-            || defined(NPY_CPU_HPPA)    \

-            || defined(NPY_CPU_PPC64)   \

-            || defined(NPY_CPU_ARMEB)   \

-            || defined(NPY_CPU_SH_BE)   \

-            || defined(NPY_CPU_MIPSEB)  \

-            || defined(NPY_CPU_OR1K)    \

-            || defined(NPY_CPU_M68K)    \

+    #elif defined(NPY_CPU_PPC)                \

+            || defined(NPY_CPU_SPARC)         \

+            || defined(NPY_CPU_S390)          \

+            || defined(NPY_CPU_HPPA)          \

+            || defined(NPY_CPU_PPC64)         \

+            || defined(NPY_CPU_ARMEB)         \

+            || defined(NPY_CPU_ARMEB_AARCH32) \

+            || defined(NPY_CPU_ARMEB_AARCH64) \

+            || defined(NPY_CPU_SH_BE)         \

+            || defined(NPY_CPU_MIPSEB)        \

+            || defined(NPY_CPU_OR1K)          \

+            || defined(NPY_CPU_M68K)          \

             || defined(NPY_CPU_ARCEB)

         #define NPY_BYTE_ORDER NPY_BIG_ENDIAN

     #else

diff --git a/numpy/core/setup.py b/numpy/core/setup.py

index 371df5b..cc1bb5b 100644

--- a/numpy/core/setup.py

+++ b/numpy/core/setup.py

@@ -452,17 +452,8 @@ def configuration(parent_package='',top_path=None):

                 moredefs.append(('NPY_RELAXED_STRIDES_DEBUG', 1))

             # Get long double representation

-            if sys.platform != 'darwin':

-                rep = check_long_double_representation(config_cmd)

-                if rep in ['INTEL_EXTENDED_12_BYTES_LE',

-                           'INTEL_EXTENDED_16_BYTES_LE',

-                           'MOTOROLA_EXTENDED_12_BYTES_BE',

-                           'IEEE_QUAD_LE', 'IEEE_QUAD_BE',

-                           'IEEE_DOUBLE_LE', 'IEEE_DOUBLE_BE',

-                           'DOUBLE_DOUBLE_BE', 'DOUBLE_DOUBLE_LE']:

-                    moredefs.append(('HAVE_LDOUBLE_%s' % rep, 1))

-                else:

-                    raise ValueError("Unrecognized long double format: %s" % rep)

+            rep = check_long_double_representation(config_cmd)

+            moredefs.append(('HAVE_LDOUBLE_%s' % rep, 1))

             # Py3K check

             if sys.version_info[0] == 3:

diff --git a/numpy/core/setup_common.py b/numpy/core/setup_common.py

index bd093c5..343318c 100644

--- a/numpy/core/setup_common.py

+++ b/numpy/core/setup_common.py

@@ -333,9 +333,9 @@ _MOTOROLA_EXTENDED_12B = ['300', '031', '000', '000', '353', '171',

 _IEEE_QUAD_PREC_BE = ['300', '031', '326', '363', '105', '100', '000', '000',

                       '000', '000', '000', '000', '000', '000', '000', '000']

 _IEEE_QUAD_PREC_LE = _IEEE_QUAD_PREC_BE[::-1]

-_DOUBLE_DOUBLE_BE = (['301', '235', '157', '064', '124', '000', '000', '000'] +

+_IBM_DOUBLE_DOUBLE_BE = (['301', '235', '157', '064', '124', '000', '000', '000'] +

                      ['000'] * 8)

-_DOUBLE_DOUBLE_LE = (['000', '000', '000', '124', '064', '157', '235', '301'] +

+_IBM_DOUBLE_DOUBLE_LE = (['000', '000', '000', '124', '064', '157', '235', '301'] +

                      ['000'] * 8)

 def long_double_representation(lines):

@@ -362,11 +362,16 @@ def long_double_representation(lines):

             # the long double

             if read[-8:] == _AFTER_SEQ:

                 saw = copy.copy(read)

+                # if the content was 12 bytes, we only have 32 - 8 - 12 = 12

+                # "before" bytes. In other words the first 4 "before" bytes went

+                # past the sliding window.

                 if read[:12] == _BEFORE_SEQ[4:]:

                     if read[12:-8] == _INTEL_EXTENDED_12B:

                         return 'INTEL_EXTENDED_12_BYTES_LE'

                     if read[12:-8] == _MOTOROLA_EXTENDED_12B:

                         return 'MOTOROLA_EXTENDED_12_BYTES_BE'

+                # if the content was 16 bytes, we are left with 32-8-16 = 16

+                # "before" bytes, so 8 went past the sliding window.

                 elif read[:8] == _BEFORE_SEQ[8:]:

                     if read[8:-8] == _INTEL_EXTENDED_16B:

                         return 'INTEL_EXTENDED_16_BYTES_LE'

@@ -374,10 +379,11 @@ def long_double_representation(lines):

                         return 'IEEE_QUAD_BE'

                     elif read[8:-8] == _IEEE_QUAD_PREC_LE:

                         return 'IEEE_QUAD_LE'

-                    elif read[8:-8] == _DOUBLE_DOUBLE_BE:

-                        return 'DOUBLE_DOUBLE_BE'

-                    elif read[8:-8] == _DOUBLE_DOUBLE_LE:

-                        return 'DOUBLE_DOUBLE_LE'

+                    elif read[8:-8] == _IBM_DOUBLE_DOUBLE_LE:

+                        return 'IBM_DOUBLE_DOUBLE_LE'

+                    elif read[8:-8] == _IBM_DOUBLE_DOUBLE_BE:

+                        return 'IBM_DOUBLE_DOUBLE_BE'

+                # if the content was 8 bytes, left with 32-8-8 = 16 bytes

                 elif read[:16] == _BEFORE_SEQ:

                     if read[16:-8] == _IEEE_DOUBLE_LE:

                         return 'IEEE_DOUBLE_LE'

diff --git a/numpy/core/src/multiarray/dragon4.c b/numpy/core/src/multiarray/dragon4.c

index e005234..e3da79f 100644

--- a/numpy/core/src/multiarray/dragon4.c

+++ b/numpy/core/src/multiarray/dragon4.c

@@ -42,6 +42,18 @@

 #define DEBUG_ASSERT(stmnt) do {} while(0)

 #endif

+static inline npy_uint64

+bitmask_u64(npy_uint32 n)

+{

+    return ~(~((npy_uint64)0) << n);

+}

+

+static inline npy_uint32

+bitmask_u32(npy_uint32 n)

+{

+    return ~(~((npy_uint32)0) << n);

+}

+

 /*

  *  Get the log base 2 of a 32-bit unsigned integer.

  *  http://graphics.stanford.edu/~seander/bithacks.html#IntegerLogLookup

@@ -102,6 +114,17 @@ LogBase2_64(npy_uint64 val)

     return LogBase2_32((npy_uint32)val);

 }

+#if defined(HAVE_LDOUBLE_IEEE_QUAD_LE) || defined(HAVE_LDOUBLE_IEEE_QUAD_BE)

+static npy_uint32

+LogBase2_128(npy_uint64 hi, npy_uint64 lo)

+{

+    if (hi) {

+        return 64 + LogBase2_64(hi);

+    }

+

+    return LogBase2_64(lo);

+}

+#endif /* HAVE_LDOUBLE_IEEE_QUAD_LE */

 /*

  * Maximum number of 32 bit blocks needed in high precision arithmetic to print

@@ -122,6 +145,45 @@ typedef struct BigInt {

     npy_uint32 blocks[c_BigInt_MaxBlocks];

 } BigInt;

+/*

+ * Dummy implementation of a memory manager for BigInts. Currently, only

+ * supports a single call to Dragon4, but that is OK because Dragon4

+ * does not release the GIL.

+ *

+ * We try to raise an error anyway if dragon4 re-enters, and this code serves

+ * as a placeholder if we want to make it re-entrant in the future.

+ *

+ * Each call to dragon4 uses 7 BigInts.

+ */

+#define BIGINT_DRAGON4_GROUPSIZE 7

+typedef struct {

+    BigInt bigints[BIGINT_DRAGON4_GROUPSIZE];

+    char repr[16384];

+} Dragon4_Scratch;

+

+static int _bigint_static_in_use = 0;

+static Dragon4_Scratch _bigint_static;

+

+static Dragon4_Scratch*

+get_dragon4_bigint_scratch(void) {

+    /* this test+set is not threadsafe, but no matter because we have GIL */

+    if (_bigint_static_in_use) {

+        PyErr_SetString(PyExc_RuntimeError,

+            "numpy float printing code is not re-entrant. "

+            "Ping the devs to fix it.");

+        return NULL;

+    }

+    _bigint_static_in_use = 1;

+

+    /* in this dummy implementation we only return the static allocation */

+    return &_bigint_static;

+}

+

+static void

+free_dragon4_bigint_scratch(Dragon4_Scratch *mem){

+    _bigint_static_in_use = 0;

+}

+

 /* Copy integer */

 static void

 BigInt_Copy(BigInt *dst, const BigInt *src)

@@ -139,32 +201,86 @@ BigInt_Copy(BigInt *dst, const BigInt *src)

 static void

 BigInt_Set_uint64(BigInt *i, npy_uint64 val)

 {

-    if (val > 0xFFFFFFFF) {

-        i->blocks[0] = val & 0xFFFFFFFF;

-        i->blocks[1] = (val >> 32) & 0xFFFFFFFF;

+    if (val > bitmask_u64(32)) {

+        i->blocks[0] = val & bitmask_u64(32);

+        i->blocks[1] = (val >> 32) & bitmask_u64(32);

         i->length = 2;

     }

     else if (val != 0) {

-        i->blocks[0] = val & 0xFFFFFFFF;

+        i->blocks[0] = val & bitmask_u64(32);

+        i->length = 1;

+    }

+    else {

+        i->length = 0;

+    }

+}

+

+#if (defined(HAVE_LDOUBLE_IBM_DOUBLE_DOUBLE_LE) || \

+     defined(HAVE_LDOUBLE_IBM_DOUBLE_DOUBLE_BE) || \

+     defined(HAVE_LDOUBLE_IEEE_QUAD_LE) || \

+     defined(HAVE_LDOUBLE_IEEE_QUAD_BE))

+static void

+BigInt_Set_2x_uint64(BigInt *i, npy_uint64 hi, npy_uint64 lo)

+{

+    if (hi > bitmask_u64(32)) {

+        i->length = 4;

+    }

+    else if (hi != 0) {

+        i->length = 3;

+    }

+    else if (lo > bitmask_u64(32)) {

+        i->length = 2;

+    }

+    else if (lo != 0) {

         i->length = 1;

     }

     else {

         i->length = 0;

     }

+

+    switch (i->length) {

+        case 4:

+            i->blocks[3] = (hi >> 32) & bitmask_u64(32);

+        case 3:

+            i->blocks[2] = hi & bitmask_u64(32);

+        case 2:

+            i->blocks[1] = (lo >> 32) & bitmask_u64(32);

+        case 1:

+            i->blocks[0] = lo & bitmask_u64(32);

+    }

 }

+#endif /* DOUBLE_DOUBLE and QUAD */

 static void

 BigInt_Set_uint32(BigInt *i, npy_uint32 val)

 {

     if (val != 0) {

         i->blocks[0] = val;

-        i->length = (val != 0);

+        i->length = 1;

     }

     else {

         i->length = 0;

     }

 }

+/*

+ * Returns 1 if the value is zero

+ */

+static int

+BigInt_IsZero(const BigInt *i)

+{

+    return i->length == 0;

+}

+

+/*

+ * Returns 1 if the value is even

+ */

+static int

+BigInt_IsEven(const BigInt *i)

+{

+    return (i->length == 0) || ( (i->blocks[0] % 2) == 0);

+}

+

 /*

  * Returns 0 if (lhs = rhs), negative if (lhs < rhs), positive if (lhs > rhs)

  */

@@ -228,7 +344,7 @@ BigInt_Add(BigInt *result, const BigInt *lhs, const BigInt *rhs)

         npy_uint64 sum = carry + (npy_uint64)(*largeCur) +

                                  (npy_uint64)(*smallCur);

         carry = sum >> 32;

-        *resultCur = sum & 0xFFFFFFFF;

+        *resultCur = sum & bitmask_u64(32);

         ++largeCur;

         ++smallCur;

         ++resultCur;

@@ -238,7 +354,7 @@ BigInt_Add(BigInt *result, const BigInt *lhs, const BigInt *rhs)

     while (largeCur != largeEnd) {

         npy_uint64 sum = carry + (npy_uint64)(*largeCur);

         carry = sum >> 32;

-        (*resultCur) = sum & 0xFFFFFFFF;

+        (*resultCur) = sum & bitmask_u64(32);

         ++largeCur;

         ++resultCur;

     }

@@ -307,13 +423,13 @@ BigInt_Multiply(BigInt *result, const BigInt *lhs, const BigInt *rhs)

                 npy_uint64 product = (*resultCur) +

                                      (*largeCur)*(npy_uint64)multiplier + carry;

                 carry = product >> 32;

-                *resultCur = product & 0xFFFFFFFF;

+                *resultCur = product & bitmask_u64(32);

                 ++largeCur;

                 ++resultCur;

             } while(largeCur != large->blocks + large->length);

             DEBUG_ASSERT(resultCur < result->blocks + maxResultLen);

-            *resultCur = (npy_uint32)(carry & 0xFFFFFFFF);

+            *resultCur = (npy_uint32)(carry & bitmask_u64(32));

         }

     }

@@ -337,7 +453,7 @@ BigInt_Multiply_int(BigInt *result, const BigInt *lhs, npy_uint32 rhs)

     const npy_uint32 *pLhsEnd = lhs->blocks + lhs->length;

     for ( ; pLhsCur != pLhsEnd; ++pLhsCur, ++resultCur) {

         npy_uint64 product = (npy_uint64)(*pLhsCur) * rhs + carry;

-        *resultCur = (npy_uint32)(product & 0xFFFFFFFF);

+        *resultCur = (npy_uint32)(product & bitmask_u64(32));

         carry = product >> 32;

     }

@@ -414,7 +530,7 @@ BigInt_Multiply10(BigInt *result)

     npy_uint32 *end = result->blocks + result->length;

     for ( ; cur != end; ++cur) {

         npy_uint64 product = (npy_uint64)(*cur) * 10ull + carry;

-        (*cur) = (npy_uint32)(product & 0xFFFFFFFF);

+        (*cur) = (npy_uint32)(product & bitmask_u64(32));

         carry = product >> 32;

     }

@@ -637,13 +753,11 @@ static BigInt g_PowerOf10_Big[] =

 /* result = 10^exponent */

 static void

-BigInt_Pow10(BigInt *result, npy_uint32 exponent)

+BigInt_Pow10(BigInt *result, npy_uint32 exponent, BigInt *temp)

 {

-    /* create two temporary values to reduce large integer copy operations */

-    BigInt temp1;

-    BigInt temp2;

-    BigInt *curTemp = &temp1;

-    BigInt *pNextTemp = &temp2;

+    /* use two temporary values to reduce large integer copy operations */

+    BigInt *curTemp = result;

+    BigInt *pNextTemp = temp;

     npy_uint32 smallExponent;

     npy_uint32 tableIdx = 0;

@@ -654,7 +768,7 @@ BigInt_Pow10(BigInt *result, npy_uint32 exponent)

      * initialize the result by looking up a 32-bit power of 10 corresponding to

      * the first 3 bits

      */

-    smallExponent = exponent & 0x7;

+    smallExponent = exponent & bitmask_u32(3);

     BigInt_Set_uint32(curTemp, g_PowerOf10_U32[smallExponent]);

     /* remove the low bits that we used for the 32-bit lookup table */

@@ -681,19 +795,17 @@ BigInt_Pow10(BigInt *result, npy_uint32 exponent)

     }

     /* output the result */

-    BigInt_Copy(result, curTemp);

+    if (curTemp != result) {

+        BigInt_Copy(result, curTemp);

+    }

 }

-/* result = in * 10^exponent */

+/* in = in * 10^exponent */

 static void

-BigInt_MultiplyPow10(BigInt *result, const BigInt *in, npy_uint32 exponent)

+BigInt_MultiplyPow10(BigInt *in, npy_uint32 exponent, BigInt *temp)

 {

-

-    /* create two temporary values to reduce large integer copy operations */

-    BigInt temp1;

-    BigInt temp2;

-    BigInt *curTemp = &temp1;

-    BigInt *pNextTemp = &temp2;

+    /* use two temporary values to reduce large integer copy operations */

+    BigInt *curTemp, *pNextTemp;

     npy_uint32 smallExponent;

     npy_uint32 tableIdx = 0;

@@ -704,12 +816,15 @@ BigInt_MultiplyPow10(BigInt *result, const BigInt *in, npy_uint32 exponent)

      * initialize the result by looking up a 32-bit power of 10 corresponding to

      * the first 3 bits

      */

-    smallExponent = exponent & 0x7;

+    smallExponent = exponent & bitmask_u32(3);

     if (smallExponent != 0) {

-        BigInt_Multiply_int(curTemp, in, g_PowerOf10_U32[smallExponent]);

+        BigInt_Multiply_int(temp, in, g_PowerOf10_U32[smallExponent]);

+        curTemp = temp;

+        pNextTemp = in;

     }

     else {

-        BigInt_Copy(curTemp, in);

+        curTemp = in;

+        pNextTemp = temp;

     }

     /* remove the low bits that we used for the 32-bit lookup table */

@@ -724,7 +839,7 @@ BigInt_MultiplyPow10(BigInt *result, const BigInt *in, npy_uint32 exponent)

             /* multiply into the next temporary */

             BigInt_Multiply(pNextTemp, curTemp, &g_PowerOf10_Big[tableIdx]);

-            // swap to the next temporary

+            /* swap to the next temporary */

             pSwap = curTemp;

             curTemp = pNextTemp;

             pNextTemp = pSwap;

@@ -736,7 +851,9 @@ BigInt_MultiplyPow10(BigInt *result, const BigInt *in, npy_uint32 exponent)

     }

     /* output the result */

-    BigInt_Copy(result, curTemp);

+    if (curTemp != in){

+        BigInt_Copy(in, curTemp);

+    }

 }

 /* result = 2^exponent */

@@ -788,7 +905,7 @@ BigInt_DivideWithRemainder_MaxQuotient9(BigInt *dividend, const BigInt *divisor)

      */

     DEBUG_ASSERT(!divisor->length == 0 &&

                 divisor->blocks[divisor->length-1] >= 8 &&

-                divisor->blocks[divisor->length-1] < 0xFFFFFFFF &&

+                divisor->blocks[divisor->length-1] < bitmask_u64(32) &&

                 dividend->length <= divisor->length);

     /*

@@ -825,10 +942,10 @@ BigInt_DivideWithRemainder_MaxQuotient9(BigInt *dividend, const BigInt *divisor)

             carry = product >> 32;

             difference = (npy_uint64)*dividendCur

-                       - (product & 0xFFFFFFFF) - borrow;

+                       - (product & bitmask_u64(32)) - borrow;

             borrow = (difference >> 32) & 1;

-            *dividendCur = difference & 0xFFFFFFFF;

+            *dividendCur = difference & bitmask_u64(32);

             ++divisorCur;

             ++dividendCur;

@@ -860,7 +977,7 @@ BigInt_DivideWithRemainder_MaxQuotient9(BigInt *dividend, const BigInt *divisor)

                                   - (npy_uint64)*divisorCur - borrow;

             borrow = (difference >> 32) & 1;

-            *dividendCur = difference & 0xFFFFFFFF;

+            *dividendCur = difference & bitmask_u64(32);

             ++divisorCur;

             ++dividendCur;

@@ -993,12 +1110,20 @@ BigInt_ShiftLeft(BigInt *result, npy_uint32 shift)

  * There is some more documentation of these changes on Ryan Juckett's website

  * at http://www.ryanjuckett.com/programming/printing-floating-point-numbers/

  *

- * Ryan Juckett's implementation did not implement "IEEE unbiased rounding",

- * except in the last digit. This has been added back, following the Burger &

- * Dybvig code, using the isEven variable.

+ * This code also has a few implementation differences from Ryan Juckett's

+ * version:

+ *  1. fixed overflow problems when mantissa was 64 bits (in float128 types),

+ *     by replacing multiplication by 2 or 4 by BigInt_ShiftLeft calls.

+ *  2. Increased c_BigInt_MaxBlocks, for 128-bit floats

+ *  3. Added more entries to the g_PowerOf10_Big table, for 128-bit floats.

+ *  4. Added unbiased rounding calculation with isEven. Ryan Juckett's

+ *     implementation did not implement "IEEE unbiased rounding", except in the

+ *     last digit. This has been added back, following the Burger & Dybvig

+ *     code, using the isEven variable.

  *

  * Arguments:

- *   * mantissa - value significand

+ *   * bigints - memory to store all bigints needed (7) for dragon4 computation.

+ *               The first BigInt should be filled in with the mantissa.

  *   * exponent - value exponent in base 2

  *   * mantissaBit - index of the highest set mantissa bit

  *   * hasUnequalMargins - is the high margin twice as large as the low margin

@@ -1007,9 +1132,11 @@ BigInt_ShiftLeft(BigInt *result, npy_uint32 shift)

  *   * pOutBuffer - buffer to output into

  *   * bufferSize - maximum characters that can be printed to pOutBuffer

  *   * pOutExponent - the base 10 exponent of the first digit

+ *

+ * Returns the number of digits written to the output buffer.

  */

 static npy_uint32

-Dragon4(const npy_uint64 mantissa, const npy_int32 exponent,

+Dragon4(BigInt *bigints, const npy_int32 exponent,

         const npy_uint32 mantissaBit, const npy_bool hasUnequalMargins,

         const DigitMode digitMode, const CutoffMode cutoffMode,

         npy_int32 cutoffNumber, char *pOutBuffer,

@@ -1025,21 +1152,24 @@ Dragon4(const npy_uint64 mantissa, const npy_int32 exponent,

      * Here, marginLow and marginHigh represent 1/2 of the distance to the next

      * floating point value above/below the mantissa.

      *

-     * scaledMarginHigh is a pointer so that it can point to scaledMarginLow in

-     * the case they must be equal to each other, otherwise it will point to

-     * optionalMarginHigh.

+     * scaledMarginHigh will point to scaledMarginLow in the case they must be

+     * equal to each other, otherwise it will point to optionalMarginHigh.

      */

-    BigInt scale;

-    BigInt scaledValue;

-    BigInt scaledMarginLow;

+    BigInt *mantissa = &bigints[0];  /* the only initialized bigint */

+    BigInt *scale = &bigints[1];

+    BigInt *scaledValue = &bigints[2];

+    BigInt *scaledMarginLow = &bigints[3];

     BigInt *scaledMarginHigh;

-    BigInt optionalMarginHigh;

+    BigInt *optionalMarginHigh = &bigints[4];

+

+    BigInt *temp1 = &bigints[5];

+    BigInt *temp2 = &bigints[6];

     const npy_float64 log10_2 = 0.30102999566398119521373889472449;

     npy_int32 digitExponent, cutoffExponent, hiBlock;

     npy_uint32 outputDigit;    /* current digit being output */

     npy_uint32 outputLen;

-    npy_bool isEven = (mantissa % 2) == 0;

+    npy_bool isEven = BigInt_IsEven(mantissa);

     npy_int32 cmp;

     /* values used to determine how to round */

@@ -1048,12 +1178,14 @@ Dragon4(const npy_uint64 mantissa, const npy_int32 exponent,

     DEBUG_ASSERT(bufferSize > 0);

     /* if the mantissa is zero, the value is zero regardless of the exponent */

-    if (mantissa == 0) {

+    if (BigInt_IsZero(mantissa)) {

         *curDigit = '0';

         *pOutExponent = 0;

         return 1;

     }

+    BigInt_Copy(scaledValue, mantissa);

+

     if (hasUnequalMargins) {

         /* if we have no fractional component */

         if (exponent > 0) {

@@ -1067,17 +1199,13 @@ Dragon4(const npy_uint64 mantissa, const npy_int32 exponent,

              */

             /* scaledValue      = 2 * 2 * mantissa*2^exponent */

-            BigInt_Set_uint64(&scaledValue, mantissa);

-            BigInt_ShiftLeft(&scaledValue, exponent + 2);

-

+            BigInt_ShiftLeft(scaledValue, exponent + 2);

             /* scale            = 2 * 2 * 1 */

-            BigInt_Set_uint32(&scale,  4);

-

+            BigInt_Set_uint32(scale,  4);

             /* scaledMarginLow  = 2 * 2^(exponent-1) */

-            BigInt_Pow2(&scaledMarginLow, exponent);

-

+            BigInt_Pow2(scaledMarginLow, exponent);

             /* scaledMarginHigh = 2 * 2 * 2^(exponent-1) */

-            BigInt_Pow2(&optionalMarginHigh, exponent + 1);

+            BigInt_Pow2(optionalMarginHigh, exponent + 1);

         }

         /* else we have a fractional exponent */

         else {

@@ -1087,34 +1215,27 @@ Dragon4(const npy_uint64 mantissa, const npy_int32 exponent,

              */

             /* scaledValue      = 2 * 2 * mantissa */

-            BigInt_Set_uint64(&scaledValue, mantissa);

-            BigInt_ShiftLeft(&scaledValue, 2);

-

+            BigInt_ShiftLeft(scaledValue, 2);

             /* scale            = 2 * 2 * 2^(-exponent) */

-            BigInt_Pow2(&scale, -exponent + 2);

-

+            BigInt_Pow2(scale, -exponent + 2);

             /* scaledMarginLow  = 2 * 2^(-1) */

-            BigInt_Set_uint32(&scaledMarginLow, 1);

-

+            BigInt_Set_uint32(scaledMarginLow, 1);

             /* scaledMarginHigh = 2 * 2 * 2^(-1) */

-            BigInt_Set_uint32(&optionalMarginHigh, 2);

+            BigInt_Set_uint32(optionalMarginHigh, 2);

         }

         /* the high and low margins are different */

-        scaledMarginHigh = &optionalMarginHigh;

+        scaledMarginHigh = optionalMarginHigh;

     }

     else {

         /* if we have no fractional component */