* mpfr/*: Update to mpfr cvs 2003-08-30.

24 files changed, 396 insertions, 154 deletions

diff --git a/mpfr/TODO b/mpfr/TODO
index ace86ff27..9f1db92c5 100644
--- a/mpfr/TODO
+++ b/mpfr/TODO
@@ -99,6 +99,8 @@ New functions to implement:
 - mpfr_set_float (MPF -> float) or mpfr_strtof (idem for mpfr_get_*)
 - implement accurate summation algorithms from Demmel
 	(http://www.cs.berkeley.edu/~demmel/AccurateSummation.ps)
+- mpfr_root to compute x^(1/n) for n integer [suggested by Damien Fisher,
+	damien@maths.usyd.edu.au, 20 Jul 2003]
 
 Rounding:
 
@@ -120,6 +122,14 @@ Efficiency:
       ...
     else
       plain operands
+
+- MPFR_CLEAR_NAN is inefficient, MPFR_SET_INF and friends should do a
+  single mask and set of their desired flags.
+
+- Nan, inf and the sign might be better off in the _mpfr_exp field,
+  since that would mean only one field set by normal operations, and
+  set by a simple store since the whole field would be updated.
+
 - mpf_t uses a scheme where the number of limbs actually present can
   be less than the selected precision, thereby allowing low precision
   values (for instance small integers) to be stored and manipulated in
@@ -172,6 +182,9 @@ Portability:
   this operation on an integer, which could be done with
   count_leading_zeros from longlong.h.
 
+- mpfr_get_d3 assumes a double is IEEE format.
+
+
 Possible future MPF / MPFR integration:
 
 - mpf routines can become "extern inline"s calling mpfr equivalents,
@@ -194,3 +207,7 @@ Possible future MPF / MPFR integration:
 
 - mpfr routines replacing mpf routines must be reentrant and thread
   safe, since of course that's what has been documented for mpf.
+
+- mpfr_random will not be wanted since there's no corresponding
+  mpf_random and new routines should not use the old style global
+  random state.
diff --git a/mpfr/acinclude.m4 b/mpfr/acinclude.m4
index eae89cf32..36aeed567 100644
--- a/mpfr/acinclude.m4
+++ b/mpfr/acinclude.m4
@@ -81,7 +81,7 @@ AC_REQUIRE([AC_HEADER_TIME])
 
 # CPU-dependent objects for the test programs
 case $host in
-  X86_PATTERN)
+  X86_PATTERN | amd64-*-*)
     AC_SUBST(TESTS_ASM_OBJECTS, x86.$OBJEXT)
     AC_DEFINE(MPFR_HAVE_TESTS_x86, 1,
               [Define to 1 if mpfr x86 test routines are available.])
diff --git a/mpfr/agm.c b/mpfr/agm.c
index 9bfde443c..aeab2702e 100644
--- a/mpfr/agm.c
+++ b/mpfr/agm.c
@@ -117,14 +117,19 @@ mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode)
   while (go_on) {
     int err, can_round;
     mp_prec_t eq;
+    double erraux;
     
-    err=1 + (int) ((3.0/2.0*(double)__gmpfr_ceil_log2((double)p)+1.0)*__gmpfr_ceil_exp2(-(double)p)
-	     +3.0*__gmpfr_ceil_exp2(-2.0*(double)p*uo/(vo-uo)));
+    erraux = (vo == uo) ? 0.0 : __gmpfr_ceil_exp2 (-2.0 * (double) p * uo 
+                                                   / (vo - uo));
+    err = 1 + (int) ((3.0 / 2.0 * (double) __gmpfr_ceil_log2 ((double) p) 
+                      + 1.0) * __gmpfr_ceil_exp2 (- (double) p)
+                     + 3.0 * erraux);
+
     if(p-err-3<=q) {
       p=q+err+4;
       err= 1 + 
 	(int) ((3.0/2.0*__gmpfr_ceil_log2((double)p)+1.0)*__gmpfr_ceil_exp2(-(double)p)
-	       +3.0*__gmpfr_ceil_exp2(-2.0*(double)p*uo/(vo-uo)));
+	       +3.0 * erraux);
     }
 
     /* Calculus of un and vn */
diff --git a/mpfr/cmp_d.c b/mpfr/cmp_d.c
index 046686605..9a0eb67c7 100644
--- a/mpfr/cmp_d.c
+++ b/mpfr/cmp_d.c
@@ -32,7 +32,7 @@ mpfr_cmp_d (mpfr_srcptr b, double d)
 
   MPFR_ASSERTN(!MPFR_IS_NAN(b));
 
-  mpfr_init2 (tmp, 53); 
+  mpfr_init2 (tmp, IEEE_DBL_MANT_DIG);
   mpfr_set_d (tmp, d, GMP_RNDN); 
 
   z = mpfr_cmp (b, tmp); 
diff --git a/mpfr/get_d.c b/mpfr/get_d.c
index 01ab28865..42526e46e 100644
--- a/mpfr/get_d.c
+++ b/mpfr/get_d.c
@@ -31,9 +31,6 @@ MA 02111-1307, USA. */
 
 static double mpfr_scale2 _PROTO ((double, int));
 
-#define IEEE_DBL_MANT_DIG 53
-
-
 /* "double" NaN and infinities are written as explicit bytes to be sure of
    getting what we want, and to be sure of not depending on libm.
 
@@ -51,24 +48,24 @@ struct dbl_bytes {
   double d;
 };
 
-#define MPFR_DBL_INFP  (* (double *) dbl_infp.b)
-#define MPFR_DBL_INFM  (* (double *) dbl_infm.b)
-#define MPFR_DBL_NAN   (* (double *) dbl_nan.b)
+#define MPFR_DBL_INFP  (* (const double *) dbl_infp.b)
+#define MPFR_DBL_INFM  (* (const double *) dbl_infm.b)
+#define MPFR_DBL_NAN   (* (const double *) dbl_nan.b)
 
 #if HAVE_DOUBLE_IEEE_LITTLE_ENDIAN
-static struct dbl_bytes dbl_infp = { { 0, 0, 0, 0, 0, 0, 0xF0, 0x7F } };
-static struct dbl_bytes dbl_infm = { { 0, 0, 0, 0, 0, 0, 0xF0, 0xFF } };
-static struct dbl_bytes dbl_nan  = { { 0, 0, 0, 0, 0, 0, 0xF8, 0x7F } };
+static const struct dbl_bytes dbl_infp = { { 0, 0, 0, 0, 0, 0, 0xF0, 0x7F } };
+static const struct dbl_bytes dbl_infm = { { 0, 0, 0, 0, 0, 0, 0xF0, 0xFF } };
+static const struct dbl_bytes dbl_nan  = { { 0, 0, 0, 0, 0, 0, 0xF8, 0x7F } };
 #endif
 #if HAVE_DOUBLE_IEEE_LITTLE_SWAPPED
-static struct dbl_bytes dbl_infp = { { 0, 0, 0xF0, 0x7F, 0, 0, 0, 0 } };
-static struct dbl_bytes dbl_infm = { { 0, 0, 0xF0, 0xFF, 0, 0, 0, 0 } };
-static struct dbl_bytes dbl_nan  = { { 0, 0, 0xF8, 0x7F, 0, 0, 0, 0 } };
+static const struct dbl_bytes dbl_infp = { { 0, 0, 0xF0, 0x7F, 0, 0, 0, 0 } };
+static const struct dbl_bytes dbl_infm = { { 0, 0, 0xF0, 0xFF, 0, 0, 0, 0 } };
+static const struct dbl_bytes dbl_nan  = { { 0, 0, 0xF8, 0x7F, 0, 0, 0, 0 } };
 #endif
 #if HAVE_DOUBLE_IEEE_BIG_ENDIAN
-static struct dbl_bytes dbl_infp = { { 0x7F, 0xF0, 0, 0, 0, 0, 0, 0 } };
-static struct dbl_bytes dbl_infm = { { 0xFF, 0xF0, 0, 0, 0, 0, 0, 0 } };
-static struct dbl_bytes dbl_nan  = { { 0x7F, 0xF8, 0, 0, 0, 0, 0, 0 } };
+static const struct dbl_bytes dbl_infp = { { 0x7F, 0xF0, 0, 0, 0, 0, 0, 0 } };
+static const struct dbl_bytes dbl_infm = { { 0xFF, 0xF0, 0, 0, 0, 0, 0, 0 } };
+static const struct dbl_bytes dbl_nan  = { { 0x7F, 0xF8, 0, 0, 0, 0, 0, 0 } };
 #endif
 #endif
 
@@ -210,9 +207,8 @@ mpfr_get_d3 (mpfr_srcptr src, mp_exp_t e, mp_rnd_t rnd_mode)
         d = 1.0;
       else
         {
-          /* Warning: the rounding may still be incorrect in the rounding
-             to the nearest mode when the result is a subnormal because of
-             a double rounding (-> 53 bits -> final precision). */
+          /* The following computations are exact thanks to the previous
+             mpfr_round_raw. */
           d = (double) tp[0] / MP_BASE_AS_DOUBLE;
           for (i = 1; i <= np; i++)
             d = (d + tp[i]) / MP_BASE_AS_DOUBLE;
@@ -247,8 +243,36 @@ mpfr_get_d1 (mpfr_srcptr src)
 }
 
 double
-mpfr_get_d_2exp (mp_exp_t *exp, mpfr_srcptr src, mp_rnd_t rnd_mode)
+mpfr_get_d_2exp (long *expptr, mpfr_srcptr src, mp_rnd_t rnd_mode)
 {
-  *exp = MPFR_GET_EXP (src);
-  return mpfr_get_d3 (src, 0, rnd_mode);
-} 
+  double ret;
+  mp_exp_t exp;
+
+  ret = mpfr_get_d3 (src, 0, rnd_mode);
+
+  if (MPFR_IS_FP(src) && MPFR_NOTZERO(src))
+    {
+      exp = MPFR_GET_EXP (src);
+
+      /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */
+      if (ret == 1.0)
+        {
+          ret = 0.5;
+          exp++;
+        }
+      else if (ret == -1.0)
+        {
+          ret = -0.5;
+          exp++;
+        }
+
+      MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0)
+                    || (ret <= -0.5 && ret > -1.0));
+      MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX);
+    }
+  else
+    exp = 0;
+
+  *expptr = exp;
+  return ret;
+}
diff --git a/mpfr/get_ld.c b/mpfr/get_ld.c
index e5837f667..0cca171ff 100644
--- a/mpfr/get_ld.c
+++ b/mpfr/get_ld.c
@@ -28,26 +28,6 @@ MA 02111-1307, USA. */
 #include "mpfr.h"
 #include "mpfr-impl.h"
 
-#ifndef LDBL_MANT_DIG
-#define LDBL_MANT_DIG 113 /* works also if long double == quad */
-#endif
-
-#ifndef DBL_MANT_DIG
-#define DBL_MANT_DIG 53
-#endif
-
-
-/* Various i386 systems have been seen with incorrect LDBL constants in
-   float.h (notes in set_ld.c), so force the value we know is right for IEEE
-   extended.  */
-
-#if HAVE_LDOUBLE_IEEE_EXT_LITTLE
-#define MPFR_LDBL_MANT_DIG   64
-#else
-#define MPFR_LDBL_MANT_DIG   LDBL_MANT_DIG
-#endif
-
-
 long double
 mpfr_get_ld (mpfr_srcptr x, mp_rnd_t rnd_mode)
 {
@@ -91,7 +71,7 @@ mpfr_get_ld (mpfr_srcptr x, mp_rnd_t rnd_mode)
         }
       /* now -1021 <= e - sh = EXP(y) <= 1023 */
       r = 0.0;
-      mpfr_init2 (z, DBL_MANT_DIG);
+      mpfr_init2 (z, IEEE_DBL_MANT_DIG);
 
       do
         {
diff --git a/mpfr/mpfr-impl.h b/mpfr/mpfr-impl.h
index 4111758fa..11ed4d4ca 100644
--- a/mpfr/mpfr-impl.h
+++ b/mpfr/mpfr-impl.h
@@ -117,6 +117,10 @@ typedef unsigned long int       mp_size_unsigned_t;
 
 /* macros for doubles, based on gmp union ieee_double_extract */
 
+#ifndef IEEE_DBL_MANT_DIG
+#define IEEE_DBL_MANT_DIG 53
+#endif
+
 typedef union ieee_double_extract Ieee_double_extract;
 
 /* for x of type ieee_double_extract */
@@ -132,6 +136,23 @@ typedef union ieee_double_extract Ieee_double_extract;
 #define DBL_NEG_INF (-1.0/0.0)
 #define DBL_NAN (0.0/0.0)
 
+/* macros for long doubles */
+
+/* we only require that LDBL_MANT_DIG is a bound on the mantissa length
+   of the "long double" type */
+#ifndef LDBL_MANT_DIG
+#define LDBL_MANT_DIG 113 /* works also if long double == quad */
+#endif
+
+/* Various i386 systems have been seen with incorrect LDBL constants in
+   float.h (notes in set_ld.c), so force the value we know is right for IEEE
+   extended.  */
+
+#if HAVE_LDOUBLE_IEEE_EXT_LITTLE
+#define MPFR_LDBL_MANT_DIG   64
+#else
+#define MPFR_LDBL_MANT_DIG   LDBL_MANT_DIG
+#endif
 
 /* LONGDOUBLE_NAN_ACTION executes the code "action" if x is a NaN. */
 
diff --git a/mpfr/mpfr-test.h b/mpfr/mpfr-test.h
index 25c98b6b5..d13533891 100644
--- a/mpfr/mpfr-test.h
+++ b/mpfr/mpfr-test.h
@@ -66,3 +66,5 @@ int ulp _PROTO ((double, double));
 double dbl _PROTO ((double, int));
 double Ulp _PROTO ((double));
 int Isnan _PROTO ((double));
+void d_trace _PROTO ((const char *, double));
+void ld_trace _PROTO ((const char *, long double));
diff --git a/mpfr/mpfr.h b/mpfr/mpfr.h
index df50a25c7..d6f5063c9 100644
--- a/mpfr/mpfr.h
+++ b/mpfr/mpfr.h
@@ -161,7 +161,7 @@ int mpfr_set_q _PROTO ((mpfr_ptr, mpq_srcptr, mp_rnd_t));
 double mpfr_get_d _PROTO ((mpfr_srcptr, mp_rnd_t));
 long double mpfr_get_ld _PROTO ((mpfr_srcptr, mp_rnd_t));
 double mpfr_get_d1 _PROTO ((mpfr_srcptr));
-double mpfr_get_d_2exp _PROTO ((mp_exp_t *, mpfr_srcptr, mp_rnd_t));
+double mpfr_get_d_2exp _PROTO ((long *, mpfr_srcptr, mp_rnd_t));
 long mpfr_get_si _PROTO ((mpfr_srcptr, mp_rnd_t));
 unsigned long mpfr_get_ui _PROTO ((mpfr_srcptr, mp_rnd_t));
 int mpfr_set_f _PROTO ((mpfr_ptr, mpf_srcptr, mp_rnd_t));
diff --git a/mpfr/mpfr.texi b/mpfr/mpfr.texi
index 7148b03f5..9c179dceb 100644
--- a/mpfr/mpfr.texi
+++ b/mpfr/mpfr.texi
@@ -901,7 +901,7 @@ Convert @var{op} to a double, using the default MPFR rounding mode
 (see function @code{mpfr_set_default_rounding_mode}).
 @end deftypefun
 
-@deftypefun double mpfr_get_d_2exp (mp_exp_t *@var{exp}, mpfr_t @var{op}, mp_rnd_t @var{rnd})
+@deftypefun double mpfr_get_d_2exp (long *@var{exp}, mpfr_t @var{op}, mp_rnd_t @var{rnd})
 Find @var{d} and @var{exp} such that @m{@var{d}\times 2^{exp}, @var{d} times 2
 raised to @var{exp}}, with @math{0.5@le{}@GMPabs{@var{d}}<1} equals
 @var{op} rounded to double precision, using the default MPFR rounding mode.
@@ -1554,7 +1554,8 @@ Assuming @var{b} is an approximation of an unknown number
 E(b)-@var{err} where E(b) is the exponent of 
 @var{b}, returns 1 if one is able to round exactly @var{x} to precision 
 @var{prec} with direction @var{rnd2},
-   and 0 otherwise. This function @strong{does not modify} its arguments.
+   and 0 otherwise (including for NaN and Inf).
+This function @strong{does not modify} its arguments.
 @end deftypefun
 
 @deftypefun mp_exp_t mpfr_get_exp (mpfr_srcptr @var{x})
diff --git a/mpfr/random2.c b/mpfr/random2.c
index fe812e958..e65d2bf66 100644
--- a/mpfr/random2.c
+++ b/mpfr/random2.c
@@ -52,7 +52,7 @@ mpfr_random2 (mpfr_ptr x, mp_size_t size, mp_exp_t exp)
   if (xn > prec + 1)
     xn = prec + 1;
 
-  /* General random mantissa.  */
+  /* Generate random mantissa.  */
   mpn_random2 (xp, xn);
 
   /* Set mandatory most significant bit.  */
diff --git a/mpfr/round_prec.c b/mpfr/round_prec.c
index 2a0a191b9..9b8bf7e17 100644
--- a/mpfr/round_prec.c
+++ b/mpfr/round_prec.c
@@ -225,10 +225,12 @@ int
 mpfr_can_round (mpfr_ptr b, mp_exp_t err, mp_rnd_t rnd1,
 		mp_rnd_t rnd2, mp_prec_t prec)
 {
-  return MPFR_IS_ZERO(b) ? 0 : /* we cannot round if b=0 */
-    mpfr_can_round_raw (MPFR_MANT(b),
-			(MPFR_PREC(b) - 1)/BITS_PER_MP_LIMB + 1,
-			MPFR_SIGN(b), err, rnd1, rnd2, prec);
+  return MPFR_IS_FP(b) ?
+    (MPFR_IS_ZERO(b) ? 0 : /* we cannot round if b=0 */
+     mpfr_can_round_raw (MPFR_MANT(b),
+			 (MPFR_PREC(b) - 1)/BITS_PER_MP_LIMB + 1,
+			 MPFR_SIGN(b), err, rnd1, rnd2, prec))
+    : 0; /* cannnot round NaN or Inf */
 }
 
 int
diff --git a/mpfr/set_d.c b/mpfr/set_d.c
index 60b08413d..a98256d3c 100644
--- a/mpfr/set_d.c
+++ b/mpfr/set_d.c
@@ -197,7 +197,7 @@ mpfr_set_d (mpfr_ptr r, double d, mp_rnd_t rnd_mode)
      since PREC(r) may be different from PREC(tmp), and then both variables
      would have same precision in the mpfr_set4 call below. */
   MPFR_MANT(tmp) = tmpmant;
-  MPFR_PREC(tmp) = 53;
+  MPFR_PREC(tmp) = IEEE_DBL_MANT_DIG;
   MPFR_SIZE(tmp) = MPFR_LIMBS_PER_DOUBLE;
 
   signd = (d < 0) ? -1 : 1;
diff --git a/mpfr/set_dfl_prec.c b/mpfr/set_dfl_prec.c
index c61875793..27ae0cd1a 100644
--- a/mpfr/set_dfl_prec.c
+++ b/mpfr/set_dfl_prec.c
@@ -24,8 +24,8 @@ MA 02111-1307, USA. */
 #include "mpfr.h"
 #include "mpfr-impl.h"
 
-/* default is 53 bits */
-mp_prec_t __gmpfr_default_fp_bit_precision = 53;
+/* default is IEEE double precision, i.e. 53 bits */
+mp_prec_t __gmpfr_default_fp_bit_precision = IEEE_DBL_MANT_DIG;
 
 void
 mpfr_set_default_prec (mp_prec_t prec)
diff --git a/mpfr/set_ld.c b/mpfr/set_ld.c
index 52954588a..743f4d8e9 100644
--- a/mpfr/set_ld.c
+++ b/mpfr/set_ld.c
@@ -21,22 +21,12 @@ the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 MA 02111-1307, USA. */
 
 #include <float.h>
-#include <limits.h>  /* for CHAR_BIT */
 
 #include "gmp.h"
 #include "gmp-impl.h"
 #include "mpfr.h"
 #include "mpfr-impl.h"
 
-#ifndef DBL_MANT_DIG
-#define DBL_MANT_DIG 53
-#endif
-
-#ifndef CHAR_BIT
-#define CHAR_BIT 8
-#endif
-
-
 /* Various i386 systems have been seen with float.h LDBL constants equal to
    the DBL ones, whereas they ought to be bigger, reflecting the 10-byte
    IEEE extended format on that processor.  gcc 3.2.1 on FreeBSD and Solaris
@@ -47,20 +37,15 @@ static const struct {
   char         bytes[10];
   long double  dummy;  /* for memory alignment */
 } ldbl_max_struct = {
-  { '\xFF','\xFF','\xFF','\xFF',
-    '\xFF','\xFF','\xFF','\xFF',
-    '\xFE','\x7F' }
+  { '\377','\377','\377','\377',
+    '\377','\377','\377','\377',
+    '\376','\177' }
 };
 #define MPFR_LDBL_MAX   (* (const long double *) ldbl_max_struct.bytes)
 #else
 #define MPFR_LDBL_MAX   LDBL_MAX
 #endif
 
-/* This is an overestimate, but fine for our purposes, it only needs to be
-   enough that "t" below can hold a long double without rounding.  */
-#define MPFR_LDBL_MANT_DIG   (CHAR_BIT * sizeof (long double))
-
-
 int
 mpfr_set_ld (mpfr_ptr r, long double d, mp_rnd_t rnd_mode)
 {
@@ -85,7 +70,7 @@ mpfr_set_ld (mpfr_ptr r, long double d, mp_rnd_t rnd_mode)
     return mpfr_set_d (r, (double) d, rnd_mode);
 
   mpfr_init2 (t, MPFR_LDBL_MANT_DIG);
-  mpfr_init2 (u, DBL_MANT_DIG);
+  mpfr_init2 (u, IEEE_DBL_MANT_DIG);
   mpfr_set_ui (t, 0, GMP_RNDN);
   while (d != 0.0)
     {
diff --git a/mpfr/tests/Makefile.am b/mpfr/tests/Makefile.am
index 0a56563bc..6fa296b6a 100644
--- a/mpfr/tests/Makefile.am
+++ b/mpfr/tests/Makefile.am
@@ -32,7 +32,7 @@ libfrtests_a_DEPENDENCIES = $(TESTS_ASM_OBJECTS)
 libfrtests_a_LIBADD = $(libfrtests_a_DEPENDENCIES)
 
 if WANT_MPFR
-check_PROGRAMS = reuse tabs tadd tagm tcan_round tcmp tcmp2 tcmp_ui tdiv tdiv_ui teq texp tfrac tget_str tlog tconst_log2 tmul tmul_2exp tmul_ui tout_str tconst_pi tpow trandom tround_prec tset_d tset_f tset_q tset_si tset_str tset_z tsqrt tsqrt_ui tui_div tui_sub tswap ttrunc trint tisnan tget_d tatan tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic texp2 tfactorial tsub tasin tconst_euler tcos tsin ttan tsub_ui tset tlog1p texpm1 tlog2 tlog10 tui_pow tpow3 tadd_ui texceptions tfma thypot tacos tgamma tset_ld tcbrt tsin_cos mpf_compat mpfr_compat tzeta tcmp_d terf
+check_PROGRAMS = reuse tabs tadd tagm tcan_round tcmp tcmp2 tcmp_ui tdiv tdiv_ui teq texp tfrac tget_str tlog tconst_log2 tmul tmul_2exp tmul_ui tout_str tconst_pi tpow trandom tround_prec tset_d tset_f tset_q tset_si tset_str tset_z tsqrt tsqrt_ui tui_div tui_sub tswap ttrunc trint tisnan tget_d tget_d_2exp tatan tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic texp2 tfactorial tsub tasin tconst_euler tcos tsin ttan tsub_ui tset tlog1p texpm1 tlog2 tlog10 tui_pow tpow3 tadd_ui texceptions tfma thypot tacos tgamma tset_ld tcbrt tsin_cos mpf_compat mpfr_compat tzeta tcmp_d terf
 check_LIBRARIES = libfrtests.a
 TESTS = $(check_PROGRAMS)
 endif
diff --git a/mpfr/tests/Makefile.in b/mpfr/tests/Makefile.in
index 37118e5ed..e48f421bc 100644
--- a/mpfr/tests/Makefile.in
+++ b/mpfr/tests/Makefile.in
@@ -206,7 +206,7 @@ EXTRA_libfrtests_a_SOURCES = x86.asm
 libfrtests_a_DEPENDENCIES = $(TESTS_ASM_OBJECTS)
 libfrtests_a_LIBADD = $(libfrtests_a_DEPENDENCIES)
 
-@WANT_MPFR_TRUE@check_PROGRAMS = reuse tabs tadd tagm tcan_round tcmp tcmp2 tcmp_ui tdiv tdiv_ui teq texp tfrac tget_str tlog tconst_log2 tmul tmul_2exp tmul_ui tout_str tconst_pi tpow trandom tround_prec tset_d tset_f tset_q tset_si tset_str tset_z tsqrt tsqrt_ui tui_div tui_sub tswap ttrunc trint tisnan tget_d tatan tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic texp2 tfactorial tsub tasin tconst_euler tcos tsin ttan tsub_ui tset tlog1p texpm1 tlog2 tlog10 tui_pow tpow3 tadd_ui texceptions tfma thypot tacos tgamma tset_ld tcbrt tsin_cos mpf_compat mpfr_compat tzeta tcmp_d terf
+@WANT_MPFR_TRUE@check_PROGRAMS = reuse tabs tadd tagm tcan_round tcmp tcmp2 tcmp_ui tdiv tdiv_ui teq texp tfrac tget_str tlog tconst_log2 tmul tmul_2exp tmul_ui tout_str tconst_pi tpow trandom tround_prec tset_d tset_f tset_q tset_si tset_str tset_z tsqrt tsqrt_ui tui_div tui_sub tswap ttrunc trint tisnan tget_d tget_d_2exp tatan tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic texp2 tfactorial tsub tasin tconst_euler tcos tsin ttan tsub_ui tset tlog1p texpm1 tlog2 tlog10 tui_pow tpow3 tadd_ui texceptions tfma thypot tacos tgamma tset_ld tcbrt tsin_cos mpf_compat mpfr_compat tzeta tcmp_d terf
 @WANT_MPFR_TRUE@check_LIBRARIES = libfrtests.a
 @WANT_MPFR_TRUE@TESTS = $(check_PROGRAMS)
 EXTRA_DIST = tgeneric.c mpf_compat.h
@@ -237,19 +237,19 @@ libfrtests_a_OBJECTS = $(am_libfrtests_a_OBJECTS)
 @WANT_MPFR_TRUE@	tsqrt_ui$(EXEEXT) tui_div$(EXEEXT) \
 @WANT_MPFR_TRUE@	tui_sub$(EXEEXT) tswap$(EXEEXT) ttrunc$(EXEEXT) \
 @WANT_MPFR_TRUE@	trint$(EXEEXT) tisnan$(EXEEXT) tget_d$(EXEEXT) \
-@WANT_MPFR_TRUE@	tatan$(EXEEXT) tcosh$(EXEEXT) tsinh$(EXEEXT) \
-@WANT_MPFR_TRUE@	ttanh$(EXEEXT) tacosh$(EXEEXT) tasinh$(EXEEXT) \
-@WANT_MPFR_TRUE@	tatanh$(EXEEXT) thyperbolic$(EXEEXT) \
-@WANT_MPFR_TRUE@	texp2$(EXEEXT) tfactorial$(EXEEXT) \
-@WANT_MPFR_TRUE@	tsub$(EXEEXT) tasin$(EXEEXT) \
-@WANT_MPFR_TRUE@	tconst_euler$(EXEEXT) tcos$(EXEEXT) \
-@WANT_MPFR_TRUE@	tsin$(EXEEXT) ttan$(EXEEXT) tsub_ui$(EXEEXT) \
-@WANT_MPFR_TRUE@	tset$(EXEEXT) tlog1p$(EXEEXT) texpm1$(EXEEXT) \
-@WANT_MPFR_TRUE@	tlog2$(EXEEXT) tlog10$(EXEEXT) tui_pow$(EXEEXT) \
-@WANT_MPFR_TRUE@	tpow3$(EXEEXT) tadd_ui$(EXEEXT) \
-@WANT_MPFR_TRUE@	texceptions$(EXEEXT) tfma$(EXEEXT) \
-@WANT_MPFR_TRUE@	thypot$(EXEEXT) tacos$(EXEEXT) tgamma$(EXEEXT) \
-@WANT_MPFR_TRUE@	tset_ld$(EXEEXT) tcbrt$(EXEEXT) \
+@WANT_MPFR_TRUE@	tget_d_2exp$(EXEEXT) tatan$(EXEEXT) \
+@WANT_MPFR_TRUE@	tcosh$(EXEEXT) tsinh$(EXEEXT) ttanh$(EXEEXT) \
+@WANT_MPFR_TRUE@	tacosh$(EXEEXT) tasinh$(EXEEXT) tatanh$(EXEEXT) \
+@WANT_MPFR_TRUE@	thyperbolic$(EXEEXT) texp2$(EXEEXT) \
+@WANT_MPFR_TRUE@	tfactorial$(EXEEXT) tsub$(EXEEXT) \
+@WANT_MPFR_TRUE@	tasin$(EXEEXT) tconst_euler$(EXEEXT) \
+@WANT_MPFR_TRUE@	tcos$(EXEEXT) tsin$(EXEEXT) ttan$(EXEEXT) \
+@WANT_MPFR_TRUE@	tsub_ui$(EXEEXT) tset$(EXEEXT) tlog1p$(EXEEXT) \
+@WANT_MPFR_TRUE@	texpm1$(EXEEXT) tlog2$(EXEEXT) tlog10$(EXEEXT) \
+@WANT_MPFR_TRUE@	tui_pow$(EXEEXT) tpow3$(EXEEXT) \
+@WANT_MPFR_TRUE@	tadd_ui$(EXEEXT) texceptions$(EXEEXT) \
+@WANT_MPFR_TRUE@	tfma$(EXEEXT) thypot$(EXEEXT) tacos$(EXEEXT) \
+@WANT_MPFR_TRUE@	tgamma$(EXEEXT) tset_ld$(EXEEXT) tcbrt$(EXEEXT) \
 @WANT_MPFR_TRUE@	tsin_cos$(EXEEXT) mpf_compat$(EXEEXT) \
 @WANT_MPFR_TRUE@	mpfr_compat$(EXEEXT) tzeta$(EXEEXT) \
 @WANT_MPFR_TRUE@	tcmp_d$(EXEEXT) terf$(EXEEXT)
@@ -457,6 +457,12 @@ tget_d_LDADD = $(LDADD)
 tget_d_DEPENDENCIES = libfrtests.a ../libmpfr.a \
 	$(top_builddir)/libgmp.la
 tget_d_LDFLAGS =
+tget_d_2exp_SOURCES = tget_d_2exp.c
+tget_d_2exp_OBJECTS = tget_d_2exp$U.$(OBJEXT)
+tget_d_2exp_LDADD = $(LDADD)
+tget_d_2exp_DEPENDENCIES = libfrtests.a ../libmpfr.a \
+	$(top_builddir)/libgmp.la
+tget_d_2exp_LDFLAGS =
 tget_str_SOURCES = tget_str.c
 tget_str_OBJECTS = tget_str$U.$(OBJEXT)
 tget_str_LDADD = $(LDADD)
@@ -699,15 +705,16 @@ DIST_SOURCES = $(libfrtests_a_SOURCES) $(EXTRA_libfrtests_a_SOURCES) \
 	tcan_round.c tcbrt.c tcmp.c tcmp2.c tcmp_d.c tcmp_ui.c \
 	tconst_euler.c tconst_log2.c tconst_pi.c tcos.c tcosh.c tdiv.c \
 	tdiv_ui.c teq.c terf.c texceptions.c texp.c texp2.c texpm1.c \
-	tfactorial.c tfma.c tfrac.c tgamma.c tget_d.c tget_str.c \
-	thyperbolic.c thypot.c tisnan.c tlog.c tlog10.c tlog1p.c \
-	tlog2.c tmul.c tmul_2exp.c tmul_ui.c tout_str.c tpow.c tpow3.c \
-	trandom.c trint.c tround_prec.c tset.c tset_d.c tset_f.c \
-	tset_ld.c tset_q.c tset_si.c tset_str.c tset_z.c tsin.c \
-	tsin_cos.c tsinh.c tsqrt.c tsqrt_ui.c tsub.c tsub_ui.c tswap.c \
-	ttan.c ttanh.c ttrunc.c tui_div.c tui_pow.c tui_sub.c tzeta.c
+	tfactorial.c tfma.c tfrac.c tgamma.c tget_d.c tget_d_2exp.c \
+	tget_str.c thyperbolic.c thypot.c tisnan.c tlog.c tlog10.c \
+	tlog1p.c tlog2.c tmul.c tmul_2exp.c tmul_ui.c tout_str.c tpow.c \
+	tpow3.c trandom.c trint.c tround_prec.c tset.c tset_d.c \
+	tset_f.c tset_ld.c tset_q.c tset_si.c tset_str.c tset_z.c \
+	tsin.c tsin_cos.c tsinh.c tsqrt.c tsqrt_ui.c tsub.c tsub_ui.c \
+	tswap.c ttan.c ttanh.c ttrunc.c tui_div.c tui_pow.c tui_sub.c \
+	tzeta.c
 DIST_COMMON = Makefile.am Makefile.in
-SOURCES = $(libfrtests_a_SOURCES) $(EXTRA_libfrtests_a_SOURCES) mpf_compat.c mpfr_compat.c reuse.c tabs.c tacos.c tacosh.c tadd.c tadd_ui.c tagm.c tasin.c tasinh.c tatan.c tatanh.c tcan_round.c tcbrt.c tcmp.c tcmp2.c tcmp_d.c tcmp_ui.c tconst_euler.c tconst_log2.c tconst_pi.c tcos.c tcosh.c tdiv.c tdiv_ui.c teq.c terf.c texceptions.c texp.c texp2.c texpm1.c tfactorial.c tfma.c tfrac.c tgamma.c tget_d.c tget_str.c thyperbolic.c thypot.c tisnan.c tlog.c tlog10.c tlog1p.c tlog2.c tmul.c tmul_2exp.c tmul_ui.c tout_str.c tpow.c tpow3.c trandom.c trint.c tround_prec.c tset.c tset_d.c tset_f.c tset_ld.c tset_q.c tset_si.c tset_str.c tset_z.c tsin.c tsin_cos.c tsinh.c tsqrt.c tsqrt_ui.c tsub.c tsub_ui.c tswap.c ttan.c ttanh.c ttrunc.c tui_div.c tui_pow.c tui_sub.c tzeta.c
+SOURCES = $(libfrtests_a_SOURCES) $(EXTRA_libfrtests_a_SOURCES) mpf_compat.c mpfr_compat.c reuse.c tabs.c tacos.c tacosh.c tadd.c tadd_ui.c tagm.c tasin.c tasinh.c tatan.c tatanh.c tcan_round.c tcbrt.c tcmp.c tcmp2.c tcmp_d.c tcmp_ui.c tconst_euler.c tconst_log2.c tconst_pi.c tcos.c tcosh.c tdiv.c tdiv_ui.c teq.c terf.c texceptions.c texp.c texp2.c texpm1.c tfactorial.c tfma.c tfrac.c tgamma.c tget_d.c tget_d_2exp.c tget_str.c thyperbolic.c thypot.c tisnan.c tlog.c tlog10.c tlog1p.c tlog2.c tmul.c tmul_2exp.c tmul_ui.c tout_str.c tpow.c tpow3.c trandom.c trint.c tround_prec.c tset.c tset_d.c tset_f.c tset_ld.c tset_q.c tset_si.c tset_str.c tset_z.c tsin.c tsin_cos.c tsinh.c tsqrt.c tsqrt_ui.c tsub.c tsub_ui.c tswap.c ttan.c ttanh.c ttrunc.c tui_div.c tui_pow.c tui_sub.c tzeta.c
 
 all: all-am
 
@@ -843,6 +850,9 @@ tgamma$(EXEEXT): $(tgamma_OBJECTS) $(tgamma_DEPENDENCIES)
 tget_d$(EXEEXT): $(tget_d_OBJECTS) $(tget_d_DEPENDENCIES) 
 	@rm -f tget_d$(EXEEXT)
 	$(LINK) $(tget_d_LDFLAGS) $(tget_d_OBJECTS) $(tget_d_LDADD) $(LIBS)
+tget_d_2exp$(EXEEXT): $(tget_d_2exp_OBJECTS) $(tget_d_2exp_DEPENDENCIES) 
+	@rm -f tget_d_2exp$(EXEEXT)
+	$(LINK) $(tget_d_2exp_LDFLAGS) $(tget_d_2exp_OBJECTS) $(tget_d_2exp_LDADD) $(LIBS)
 tget_str$(EXEEXT): $(tget_str_OBJECTS) $(tget_str_DEPENDENCIES) 
 	@rm -f tget_str$(EXEEXT)
 	$(LINK) $(tget_str_LDFLAGS) $(tget_str_OBJECTS) $(tget_str_LDADD) $(LIBS)
@@ -1065,6 +1075,8 @@ tgamma_.c: tgamma.c $(ANSI2KNR)
 	$(CPP) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) `if test -f $(srcdir)/tgamma.c; then echo $(srcdir)/tgamma.c; else echo tgamma.c; fi` | sed 's/^# \([0-9]\)/#line \1/' | $(ANSI2KNR) > $@ || rm -f $@
 tget_d_.c: tget_d.c $(ANSI2KNR)
 	$(CPP) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) `if test -f $(srcdir)/tget_d.c; then echo $(srcdir)/tget_d.c; else echo tget_d.c; fi` | sed 's/^# \([0-9]\)/#line \1/' | $(ANSI2KNR) > $@ || rm -f $@
+tget_d_2exp_.c: tget_d_2exp.c $(ANSI2KNR)
+	$(CPP) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) `if test -f $(srcdir)/tget_d_2exp.c; then echo $(srcdir)/tget_d_2exp.c; else echo tget_d_2exp.c; fi` | sed 's/^# \([0-9]\)/#line \1/' | $(ANSI2KNR) > $@ || rm -f $@
 tget_str_.c: tget_str.c $(ANSI2KNR)
 	$(CPP) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) `if test -f $(srcdir)/tget_str.c; then echo $(srcdir)/tget_str.c; else echo tget_str.c; fi` | sed 's/^# \([0-9]\)/#line \1/' | $(ANSI2KNR) > $@ || rm -f $@
 thyperbolic_.c: thyperbolic.c $(ANSI2KNR)
@@ -1163,26 +1175,26 @@ tests_.$(OBJEXT) tests_.lo texceptions_.$(OBJEXT) texceptions_.lo \
 texp_.$(OBJEXT) texp_.lo texp2_.$(OBJEXT) texp2_.lo texpm1_.$(OBJEXT) \
 texpm1_.lo tfactorial_.$(OBJEXT) tfactorial_.lo tfma_.$(OBJEXT) \
 tfma_.lo tfrac_.$(OBJEXT) tfrac_.lo tgamma_.$(OBJEXT) tgamma_.lo \
-tget_d_.$(OBJEXT) tget_d_.lo tget_str_.$(OBJEXT) tget_str_.lo \
-thyperbolic_.$(OBJEXT) thyperbolic_.lo thypot_.$(OBJEXT) thypot_.lo \
-tisnan_.$(OBJEXT) tisnan_.lo tlog_.$(OBJEXT) tlog_.lo tlog10_.$(OBJEXT) \
-tlog10_.lo tlog1p_.$(OBJEXT) tlog1p_.lo tlog2_.$(OBJEXT) tlog2_.lo \
-tmul_.$(OBJEXT) tmul_.lo tmul_2exp_.$(OBJEXT) tmul_2exp_.lo \
-tmul_ui_.$(OBJEXT) tmul_ui_.lo tout_str_.$(OBJEXT) tout_str_.lo \
-tpow_.$(OBJEXT) tpow_.lo tpow3_.$(OBJEXT) tpow3_.lo trandom_.$(OBJEXT) \
-trandom_.lo trint_.$(OBJEXT) trint_.lo tround_prec_.$(OBJEXT) \
-tround_prec_.lo tset_.$(OBJEXT) tset_.lo tset_d_.$(OBJEXT) tset_d_.lo \
-tset_f_.$(OBJEXT) tset_f_.lo tset_ld_.$(OBJEXT) tset_ld_.lo \
-tset_q_.$(OBJEXT) tset_q_.lo tset_si_.$(OBJEXT) tset_si_.lo \
-tset_str_.$(OBJEXT) tset_str_.lo tset_z_.$(OBJEXT) tset_z_.lo \
-tsin_.$(OBJEXT) tsin_.lo tsin_cos_.$(OBJEXT) tsin_cos_.lo \
-tsinh_.$(OBJEXT) tsinh_.lo tsqrt_.$(OBJEXT) tsqrt_.lo \
-tsqrt_ui_.$(OBJEXT) tsqrt_ui_.lo tsub_.$(OBJEXT) tsub_.lo \
-tsub_ui_.$(OBJEXT) tsub_ui_.lo tswap_.$(OBJEXT) tswap_.lo \
-ttan_.$(OBJEXT) ttan_.lo ttanh_.$(OBJEXT) ttanh_.lo ttrunc_.$(OBJEXT) \
-ttrunc_.lo tui_div_.$(OBJEXT) tui_div_.lo tui_pow_.$(OBJEXT) \
-tui_pow_.lo tui_sub_.$(OBJEXT) tui_sub_.lo tzeta_.$(OBJEXT) tzeta_.lo : \
-$(ANSI2KNR)
+tget_d_.$(OBJEXT) tget_d_.lo tget_d_2exp_.$(OBJEXT) tget_d_2exp_.lo \
+tget_str_.$(OBJEXT) tget_str_.lo thyperbolic_.$(OBJEXT) thyperbolic_.lo \
+thypot_.$(OBJEXT) thypot_.lo tisnan_.$(OBJEXT) tisnan_.lo \
+tlog_.$(OBJEXT) tlog_.lo tlog10_.$(OBJEXT) tlog10_.lo tlog1p_.$(OBJEXT) \
+tlog1p_.lo tlog2_.$(OBJEXT) tlog2_.lo tmul_.$(OBJEXT) tmul_.lo \
+tmul_2exp_.$(OBJEXT) tmul_2exp_.lo tmul_ui_.$(OBJEXT) tmul_ui_.lo \
+tout_str_.$(OBJEXT) tout_str_.lo tpow_.$(OBJEXT) tpow_.lo \
+tpow3_.$(OBJEXT) tpow3_.lo trandom_.$(OBJEXT) trandom_.lo \
+trint_.$(OBJEXT) trint_.lo tround_prec_.$(OBJEXT) tround_prec_.lo \
+tset_.$(OBJEXT) tset_.lo tset_d_.$(OBJEXT) tset_d_.lo tset_f_.$(OBJEXT) \
+tset_f_.lo tset_ld_.$(OBJEXT) tset_ld_.lo tset_q_.$(OBJEXT) tset_q_.lo \
+tset_si_.$(OBJEXT) tset_si_.lo tset_str_.$(OBJEXT) tset_str_.lo \
+tset_z_.$(OBJEXT) tset_z_.lo tsin_.$(OBJEXT) tsin_.lo \
+tsin_cos_.$(OBJEXT) tsin_cos_.lo tsinh_.$(OBJEXT) tsinh_.lo \
+tsqrt_.$(OBJEXT) tsqrt_.lo tsqrt_ui_.$(OBJEXT) tsqrt_ui_.lo \
+tsub_.$(OBJEXT) tsub_.lo tsub_ui_.$(OBJEXT) tsub_ui_.lo \
+tswap_.$(OBJEXT) tswap_.lo ttan_.$(OBJEXT) ttan_.lo ttanh_.$(OBJEXT) \
+ttanh_.lo ttrunc_.$(OBJEXT) ttrunc_.lo tui_div_.$(OBJEXT) tui_div_.lo \
+tui_pow_.$(OBJEXT) tui_pow_.lo tui_sub_.$(OBJEXT) tui_sub_.lo \
+tzeta_.$(OBJEXT) tzeta_.lo : $(ANSI2KNR)
 
 mostlyclean-libtool:
 	-rm -f *.lo
diff --git a/mpfr/tests/tagm.c b/mpfr/tests/tagm.c
index 57f5c6a37..100317a35 100644
--- a/mpfr/tests/tagm.c
+++ b/mpfr/tests/tagm.c
@@ -82,15 +82,32 @@ check_large (void)
 {
   mpfr_t a, b, agm;
 
-  mpfr_init2(a, 82); mpfr_init2(b, 82); mpfr_init2(agm, 82);
-  mpfr_set_ui(a, 1, GMP_RNDN);
-  mpfr_set_str_raw(b, "0.1111101100001000000001011000110111101000001011111000100001000101010100011111110010E-39");
-  mpfr_agm(agm, a, b, GMP_RNDN);
-  mpfr_set_str_raw(a, "0.1110001000111101101010101010101101001010001001001011100101111011110101111001111100E-4");
-  if (mpfr_cmp(agm, a)) {
-    fprintf(stderr, "mpfr_agm failed for precision 82\n"); exit(1);
-  }
-  mpfr_clear(a); mpfr_clear(b); mpfr_clear(agm);
+  mpfr_init2 (a, 82);
+  mpfr_init2 (b, 82);
+  mpfr_init2 (agm, 82);
+
+  mpfr_set_ui (a, 1, GMP_RNDN);
+  mpfr_set_str_raw (b, "0.1111101100001000000001011000110111101000001011111000100001000101010100011111110010E-39");
+  mpfr_agm (agm, a, b, GMP_RNDN);
+  mpfr_set_str_raw (a, "0.1110001000111101101010101010101101001010001001001011100101111011110101111001111100E-4");
+  if (mpfr_cmp (agm, a))
+    {
+      fprintf (stderr, "mpfr_agm failed for precision 82\n");
+      exit (1);
+    }
+  
+  /* problem found by Damien Fischer <damien@maths.usyd.edu.au> 4 Aug 2003:
+     produced a division by zero exception */
+  mpfr_set_prec (a, 268);
+  mpfr_set_prec (b, 268);
+  mpfr_set_prec (agm, 268);
+  mpfr_set_str (a, "703.93543315330225238487276503953366664991725089988315253092140138947103394917006", 10, GMP_RNDN);
+  mpfr_set_str (b, "703.93543315330225238487279020523738740563816490895994499256063816906728642622316", 10, GMP_RNDN);
+  mpfr_agm (agm, a, b, GMP_RNDN);
+
+  mpfr_clear (a);
+  mpfr_clear (b);
+  mpfr_clear (agm);
 }
 
 #if 0
diff --git a/mpfr/tests/tests.c b/mpfr/tests/tests.c
index a0c31cd6d..78f90b414 100644
--- a/mpfr/tests/tests.c
+++ b/mpfr/tests/tests.c
@@ -126,11 +126,6 @@ void
 tests_rand_end (void)
 {
   RANDS_CLEAR ();
-  if (__gmp_rands_initialized)
-    {
-      gmp_randclear (__gmp_rands);
-      __gmp_rands_initialized = 0;
-    }
 }
 
 /* initialization function for tests using the hardware floats */
@@ -180,7 +175,7 @@ mpfr_test_init ()
    exponent range is only enforced when storing to memory.
 
    For reference, on most i386 systems the default is 64-bit "long double"
-   precision, but on FreeBSD 3.x it's 53-bit "double".  */
+   precision, but on FreeBSD 3.x and amd64 5.x it's 53-bit "double".  */
 
 void
 tests_machine_prec_double (void)
@@ -264,3 +259,49 @@ Isnan (double d)
 {
   return (d) != (d);
 }
+
+void
+d_trace (const char *name, double d)
+{
+  union {
+    double         d;
+    unsigned char  b[sizeof(double)];
+  } u;
+  int  i;
+
+  if (name != NULL && name[0] != '\0')
+    printf ("%s=", name);
+
+  u.d = d;
+  printf ("[");
+  for (i = 0; i < sizeof (u.b); i++)
+    {
+      if (i != 0)
+        printf (" ");
+      printf ("%02X", (int) u.b[i]);
+    }
+  printf ("] %.20g\n", d);
+}
+
+void
+ld_trace (const char *name, long double ld)
+{
+  union {
+    long double    ld;
+    unsigned char  b[sizeof(long double)];
+  } u;
+  int  i;
+
+  if (name != NULL && name[0] != '\0')
+    printf ("%s=", name);
+
+  u.ld = ld;
+  printf ("[");
+  for (i = 0; i < sizeof (u.b); i++)
+    {
+      if (i != 0)
+        printf (" ");
+      printf ("%02X", (int) u.b[i]);
+    }
+  printf ("] %.20Lg\n", ld);
+}
diff --git a/mpfr/tests/tget_d_2exp.c b/mpfr/tests/tget_d_2exp.c
new file mode 100644
index 000000000..2a7b42fcc
--- /dev/null
+++ b/mpfr/tests/tget_d_2exp.c
@@ -0,0 +1,123 @@
+/* Test mpfr_get_d_2exp.
+
+Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The MPFR Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the MPFR Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <float.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+#include "mpfr-impl.h"
+#include "mpfr-test.h"
+
+
+/* Check that hardware rounding doesn't make mpfr_get_d_2exp return a value
+   outside its defined range. */
+static void
+check_round (void)
+{
+  static const unsigned long data[] = { 1, 32, 53, 54, 64, 128, 256, 512 };
+  mpfr_t  f;
+  double  got;
+  long    got_exp;
+  int     i, rnd_mode, neg;
+
+  mpfr_init2 (f, 1024L);
+
+  for (rnd_mode = 0; rnd_mode < 4; rnd_mode++)
+    {
+      for (i = 0; i < numberof (data); i++)
+        {
+          mpfr_set_ui (f, 1L, GMP_RNDZ);
+          mpfr_mul_2exp (f, f, data[i], GMP_RNDZ);
+          mpfr_sub_ui (f, f, 1L, GMP_RNDZ);
+
+          for (neg = 0; neg <= 1; neg++)
+            {
+              got = mpfr_get_d_2exp (&got_exp, f, rnd_mode);
+
+              if (neg == 0
+                  ? (got < 0.5 || got >= 1.0)
+                  : (got <= -1.0 || got > -0.5))
+                {
+                  printf  ("mpfr_get_d_2exp wrong on 2**%lu-1\n", data[i]);
+                  printf  ("result out of range, expect 0.5 <= got < 1.0\n");
+                  printf  ("   rnd_mode = %d\n", rnd_mode);
+                  printf  ("   data[i]  = %lu\n", data[i]);
+                  printf  ("   f    ");
+                  mpfr_out_str (stdout, 2, 0, f, GMP_RNDN);
+                  printf  ("\n");
+                  d_trace ("   got  ", got);
+                  printf  ("   got exp  %ld\n", got_exp);
+                  abort();
+                }
+
+              mpfr_neg (f, f, GMP_RNDZ);
+            }
+        }
+    }
+
+  mpfr_clear (f);
+}
+
+
+static void
+check_inf_nan ()
+{
+  /* only if nans and infs are available */
+#if _GMP_IEEE_FLOATS
+  mpfr_t  x;
+  double  d;
+  long    exp;
+
+  mpfr_init2 (x, 123);
+
+  mpfr_set_inf (x, 1);
+  d = mpfr_get_d_2exp (&exp, x, GMP_RNDZ);
+  ASSERT_ALWAYS (d > 0);
+  ASSERT_ALWAYS (DOUBLE_ISINF (d));
+
+  mpfr_set_inf (x, -1);
+  d = mpfr_get_d_2exp (&exp, x, GMP_RNDZ);
+  ASSERT_ALWAYS (d < 0);
+  ASSERT_ALWAYS (DOUBLE_ISINF (d));
+
+  mpfr_set_nan (x);
+  d = mpfr_get_d_2exp (&exp, x, GMP_RNDZ);
+  ASSERT_ALWAYS (DOUBLE_ISNAN (d));
+
+  mpfr_clear (x);
+#endif
+}
+
+
+int
+main (void)
+{
+  tests_start_mpfr ();
+  mpfr_test_init ();
+
+  check_round ();
+  check_inf_nan ();
+
+  tests_end_mpfr ();
+  return 0;
+}
diff --git a/mpfr/tests/tget_str.c b/mpfr/tests/tget_str.c
index 4aa0ddd9a..97e3f1d56 100644
--- a/mpfr/tests/tget_str.c
+++ b/mpfr/tests/tget_str.c
@@ -188,7 +188,7 @@ static void
 check_large (void)
 {
   mpfr_t x;
-  char *s, s1[5];
+  char *s, s1[7];
   mp_exp_t e;
 
   mpfr_init2 (x, 3322);
diff --git a/mpfr/tests/tset_ld.c b/mpfr/tests/tset_ld.c
index c0b05a1d0..947723776 100644
--- a/mpfr/tests/tset_ld.c
+++ b/mpfr/tests/tset_ld.c
@@ -59,8 +59,13 @@ check_set_get (long double d, mpfr_t x)
       e = mpfr_get_ld (x, r);
       if (e != d && !(Isnan_ld(e) && Isnan_ld(d)))
         {
-          fprintf (stderr, "Error: mpfr_get_ld o mpfr_set_ld <> Id\n");
-          fprintf (stderr, "d=%1.30Le get_ld(set_ld(d))=%1.30Le\n", d, e);
+          printf ("Error: mpfr_get_ld o mpfr_set_ld <> Id\n");
+          printf ("  r=%d\n", r);
+          printf ("  d=%1.30Le get_ld(set_ld(d))=%1.30Le\n", d, e);
+          ld_trace ("  d", d);
+          printf ("  x="); mpfr_out_str (NULL, 16, 0, x, GMP_RNDN);
+          printf ("\n");
+          ld_trace ("  e", e);
           exit (1);
         }
     }
@@ -122,11 +127,15 @@ main (int argc, char *argv[])
   check_set_get (d, x);
   check_set_get (-d, x);
 
-  /* checks 2^1024 */
+  /* checks largest 2^(2^k) that is representable as a double */
   f = 1.3407807929942597100e155; /* 2^512 */
+  i = 512;
   d = (long double) f;
-  if (sizeof(long double) > sizeof(double))
-    d = d * d; /* 2^1024 */
+  while (2 * i < LDBL_MAX_EXP)
+    {
+      d = d * d;
+      i = 2 * i;
+    }
   check_set_get (d, x);
 
   /* checks that 2^i, 2^i+1 and 2^i-1 are correctly converted */
diff --git a/mpfr/ui_pow_ui.c b/mpfr/ui_pow_ui.c
index d04852704..b8a512fb4 100644
--- a/mpfr/ui_pow_ui.c
+++ b/mpfr/ui_pow_ui.c
@@ -1,6 +1,6 @@
 /*  mpfr_ui_pow_ui -- compute the power beetween two machine integer
 
-Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
+Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc.
 
 This file is part of the MPFR Library.
 
@@ -20,6 +20,7 @@ the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 MA 02111-1307, USA. */
 
 #include "gmp.h"
+#include "gmp-impl.h"
 #include "mpfr.h"
 #include "mpfr-impl.h"
 
@@ -63,7 +64,7 @@ mpfr_ui_pow_ui (mpfr_ptr x, unsigned long int y, unsigned long int n,
 	  if (mpfr_mul (res, res, res, GMP_RNDU))
 	    inexact = 1;
 	  err++;
-	  if (n & (1 << i))
+	  if (n & (1UL << i))
 	    if (mpfr_mul_ui (res, res, y, GMP_RNDU))
 	      inexact = 1;
 	}
diff --git a/mpfr/urandomb.c b/mpfr/urandomb.c
index 423359c8e..2ac85ec01 100644
--- a/mpfr/urandomb.c
+++ b/mpfr/urandomb.c
@@ -34,6 +34,7 @@ mpfr_urandomb (mpfr_ptr rop, gmp_randstate_t rstate)
   mp_ptr rp;
   mp_prec_t nbits;
   mp_size_t nlimbs;
+  mp_size_t k; /* number of high zero limbs */
   mp_exp_t exp;
   int cnt;
 
@@ -43,16 +44,19 @@ mpfr_urandomb (mpfr_ptr rop, gmp_randstate_t rstate)
   nbits = MPFR_PREC(rop);
   nlimbs = (nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
 
-  _gmp_rand (rp, rstate, nbits);
+  _gmp_rand (rp, rstate, nlimbs * BITS_PER_MP_LIMB);
 
-  /* If nbits isn't a multiple of BITS_PER_MP_LIMB, shift up.  */
-  if (nbits % BITS_PER_MP_LIMB != 0)
-    mpn_lshift (rp, rp, nlimbs, BITS_PER_MP_LIMB - nbits % BITS_PER_MP_LIMB);
+  /* If nbits isn't a multiple of BITS_PER_MP_LIMB, mask the low bits */
+  cnt = nlimbs * BITS_PER_MP_LIMB - nbits;
+  if (cnt != 0)
+    rp[0] &= ~((MP_LIMB_T_ONE << cnt) - MP_LIMB_T_ONE);
 
   exp = 0;
+  k = 0;
   while (nlimbs != 0 && rp[nlimbs - 1] == 0)
     {
-      nlimbs--;
+      k ++;
+      nlimbs --;
       exp -= BITS_PER_MP_LIMB;
     }
 
@@ -66,11 +70,9 @@ mpfr_urandomb (mpfr_ptr rop, gmp_randstate_t rstate)
           return 1;
         }
       if (cnt != 0)
-        mpn_lshift (rp, rp, nlimbs, cnt);
-
-      cnt = (mp_prec_t) nlimbs * BITS_PER_MP_LIMB - nbits;
-      /* cnt is the number of non significant bits in the low limb */
-      rp[0] &= ~((MP_LIMB_T_ONE << cnt) - 1);
+        mpn_lshift (rp + k, rp, nlimbs, cnt);
+      if (k)
+        MPN_ZERO (rp, k);
     }
 
   MPFR_SET_POS (rop);