mini-gmp and mpq

[Marco Bodrato]

Ciao,

I reply to an e-mail from gmp-discuss.

Il Lun, 19 Febbraio 2018 5:24 pm, Marc Glisse ha scritto:
> On Mon, 19 Feb 2018, Marco Bodrato wrote:
>> Il Mer, 14 Febbraio 2018 8:07 pm, Andreas Fabri ha scritto:
>>> 2) the conversion of a mpq to its closest double
>>
>> mpq_get_d currently gives a double, but rounded towards zero...
>
> Except when it uses the broken generic code :-(

Maybe things get better if we mask a fixed number of bits during conversion?

A possible patch (to mpn/generic/get_d.c) follows:

diff -r 110bcd4c29f4 mpn/generic/get_d.c
--- a/mpn/generic/get_d.c       Sun May 06 08:48:36 2018 +0200
+++ b/mpn/generic/get_d.c       Sun May 06 18:28:20 2018 +0200
@@ -32,6 +32,12 @@
 GNU Lesser General Public License along with the GNU MP Library.  If not,
 see https://www.gnu.org/licenses/.  */

+#include "config.h"
+
+#if HAVE_FLOAT_H
+#include <float.h>  /* for DBL_MANT_DIG and FLT_RADIX */
+#endif
+
 #include "gmp-impl.h"
 #include "longlong.h"

@@ -354,7 +360,20 @@
 #endif

 #if ! FORMAT_RECOGNIZED
-    {      /* Non-IEEE or strange limb size, do something generic. */
+
+#if !defined(GMP_DBL_MANT_BITS)
+#if defined(DBL_MANT_DIG) && FLT_RADIX == 2
+#define GMP_DBL_MANT_BITS DBL_MANT_DIG
+#else
+/* FIXME: Chose a smarter default value. */
+#define GMP_DBL_MANT_BITS (16 * sizeof (double))
+#endif
+#endif
+
+    { /* Non-IEEE or strange limb size, generically convert
+        GMP_DBL_MANT_BITS bits. */
+      mp_limb_t l;
+      int m;
       mp_size_t i;
       double d, weight;
       unsigned long uexp;
@@ -366,13 +385,21 @@
         limbs.  */
-      weight = 1/MP_BASE_AS_DOUBLE;
-      d = up[size - 1];
-      for (i = size - 2; i >= 0; i--)
-       {
-         d += up[i] * weight;
-         weight /= MP_BASE_AS_DOUBLE;
-         if (weight == 0)
-           break;
-       }
+      i = size - 1;
+      l = up[i];
+      count_leading_zeros (m, l);
+      m = m - GMP_LIMB_BITS + GMP_DBL_MANT_BITS;
+      if (m < 0)
+       l &= GMP_NUMB_MAX << -m;
+      d = l;
+      for (weight = 1/MP_BASE_AS_DOUBLE; m > 0 && --i >= 0;) {
+       l = up[i];
+       m -= GMP_NUMB_BITS;
+       if (m < 0)
+         l &= GMP_NUMB_MAX << -m;
+       d += l * weight;
+       weight /= MP_BASE_AS_DOUBLE;
+       if (weight == 0)
+         break;
+      }

       /* Now apply exp.  */
       exp -= GMP_NUMB_BITS;


The large default value for GMP_DBL_MANT_BITS should give us basically the
same behaviour as the current code, in the worst case.

Ĝis,
m

-- 
http://bodrato.it/papers/