* mpfr/*: Update to 20010215.

33 files changed, 357 insertions, 390 deletions

diff --git a/mpfr/NEWS b/mpfr/NEWS
index 1792ad624..e98cf6261 100644
--- a/mpfr/NEWS
+++ b/mpfr/NEWS
@@ -14,7 +14,7 @@ Changes from version 1.0 to version 2001:
 - mpfr_init and mpfr_init2 now initialize the corresponding variable to 0
   (like in other initialization functions from GNU MP)
 - in case memory allocation fails, an error message is output
-- several bugs of version 0.4 were fixed
+- several bugs of version 1.0 were fixed
 
 Changes from version 0.4 to version 1.0:
 
diff --git a/mpfr/TODO b/mpfr/TODO
index 4c62eaf4c..2694694fb 100644
--- a/mpfr/TODO
+++ b/mpfr/TODO
@@ -16,3 +16,24 @@ Miscellaneous:
 - add mpfr_get_ld for 'long double' [asked by J-C Fauge`re] ?
   (exists since K&R, but has not necessarily a greater precision than double)
   cf http://anubis.dkuug.dk/jtc1/sc22/wg14/www/docs/n869/
+
+- from Kevin Ryde <user42@zip.com.au>:
+   Also for pi.c, a pre-calculated compiled-in pi to a few thousand
+   digits would be good value I think.  After all, say 10000 bits using
+   1250 bytes would still be small compared to the code size!
+
+- from Kevin Ryde <user42@zip.com.au>:
+   Determine the exp2/exp3 thresholds using tune/tuneup.c.
+
+- add a description of the algorithms used
+
+- add undef's in mpf2mpfr.h to avoid warnings
+
+- rewrite mpfr_get_d2 to avoid rounding problems, for example on the IA64:
+x=8.94446216651049664713e+94,74 y=-8.69096029476006043478e+78,113 pz=102,rnd=GMP_RNDU
+got 8.94446216651049812927e+94
+result should be 8.94446216651049516499e+94 (diff=-2 ulp)
+i.e.0.101010111000010100110011011001011100111011100100010110000000000000000000000000000000000000000000000000[00000000000000000000000000]E316
+
+- problem when reading a float followed by a character, for example 1.5*x
+	[from Fabrice.Rouillier@loria.fr, Mon, 04 Dec 2000]
diff --git a/mpfr/add.c b/mpfr/add.c
index e24054924..6aaf88ad7 100644
--- a/mpfr/add.c
+++ b/mpfr/add.c
@@ -29,7 +29,7 @@ MA 02111-1307, USA. */
 
 extern void mpfr_sub1 _PROTO((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, 
 			      mp_rnd_t, int));
-void mpfr_add1 _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t, int));
+void mpfr_add1 _PROTO((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t, int));
 
 #define ONE ((mp_limb_t) 1)
 
diff --git a/mpfr/agm.c b/mpfr/agm.c
index 30ab45c1b..27a992769 100644
--- a/mpfr/agm.c
+++ b/mpfr/agm.c
@@ -55,11 +55,12 @@ _mpfr_ceil_exp2 (double d)
 {
   long exp;
   union ieee_double_extract x;
-  
+
   exp = (long) d;
   if (d != (double) exp) exp++;
   /* now exp = ceil(d) */
   x.d = 1.0;
+  if (exp < -1022) exp = -1022;
   x.s.exp = 1023 + exp;
   return x.d;
 }
@@ -147,12 +148,12 @@ mpfr_agm (r, op2, op1, rnd_mode)
     
     eq=0;
     
-    err=1 + (int) ((3.0/2.0*_mpfr_ceil_log2(p)+1.0)*_mpfr_ceil_exp2(-(double)p)
+    err=1 + (int) ((3.0/2.0*(double)_mpfr_ceil_log2((double)p)+1.0)*_mpfr_ceil_exp2(-(double)p)
 	     +3.0*_mpfr_ceil_exp2(-2.0*(double)p*uo/(vo-uo)));
     if(p-err-3<=q) {
       p=q+err+4;
       err= 1 + 
-	(int) ((3.0/2.0*_mpfr_ceil_log2(p)+1.0)*_mpfr_ceil_exp2(-(double)p)
+	(int) ((3.0/2.0*_mpfr_ceil_log2((double)p)+1.0)*_mpfr_ceil_exp2(-(double)p)
 	       +3.0*_mpfr_ceil_exp2(-2.0*(double)p*uo/(vo-uo)));
     }
 
diff --git a/mpfr/clear.c b/mpfr/clear.c
index 5f492badf..0b1db32fd 100644
--- a/mpfr/clear.c
+++ b/mpfr/clear.c
@@ -34,5 +34,5 @@ mpfr_clear (m)
 #endif
 {
   /* be careful to always free an entire number of limbs */
-  (*_mp_free_func) (MPFR_MANT(m), MPFR_ABSSIZE(m) * BYTES_PER_MP_LIMB);
+  (*__gmp_free_func) (MPFR_MANT(m), MPFR_ABSSIZE(m) * BYTES_PER_MP_LIMB);
 }
diff --git a/mpfr/cmp.c b/mpfr/cmp.c
index d0cb175db..5b7b73cf5 100644
--- a/mpfr/cmp.c
+++ b/mpfr/cmp.c
@@ -185,7 +185,10 @@ mpfr_cmp2(b, c)
 	  
       /* bn < 0; if some limb of c is nonzero, return k+1, otherwise return k*/
 
-      if (cn>=0 && (cp[cn--] << (BITS_PER_MP_LIMB - d))) { return k+1; }
+      /* if d < BITS_PER_MP_LIMB, only the last d bits of cp[cn] have to be
+	 considered, otherwise all bits */
+      if (d < BITS_PER_MP_LIMB)
+	if (cn>=0 && (cp[cn--] << (BITS_PER_MP_LIMB - d))) { return k+1; }
 
       while (cn >= 0) 
 	if (cp[cn--]) return k+1; 
diff --git a/mpfr/dump.c b/mpfr/dump.c
index a0319c015..a7180bcac 100644
--- a/mpfr/dump.c
+++ b/mpfr/dump.c
@@ -59,5 +59,5 @@ mpfr_dump (u, rnd_mode)
     printf ("-0.%se%ld\n", str + 1, exp);
   else
     printf ("0.%se%ld\n", str, exp);
-  (*_mp_free_func) (str, strlen(str)+1);
+  (*__gmp_free_func) (str, strlen(str)+1);
 }
diff --git a/mpfr/generic.c b/mpfr/generic.c
index cd3b257f7..88e509367 100644
--- a/mpfr/generic.c
+++ b/mpfr/generic.c
@@ -55,30 +55,19 @@ GENERIC (y, p, r, m)
   mpz_t* qtoj;
   mpfr_t tmp;
 #endif
-  int diff,expo;
+  int diff, expo;
   int precy = MPFR_PREC(y);
   MPFR_CLEAR_FLAGS(y); 
   n = 1 << m;
-  P = (mpz_t*) (*_mp_allocate_func) ((m+1) * sizeof(mpz_t));
-  S = (mpz_t*) (*_mp_allocate_func) ((m+1) * sizeof(mpz_t));
-  ptoj = (mpz_t*) (*_mp_allocate_func) ((m+1) * sizeof(mpz_t)); /* ptoj[i] = mantissa^(2^i) */
+  P = (mpz_t*) (*__gmp_allocate_func) ((m+1) * sizeof(mpz_t));
+  S = (mpz_t*) (*__gmp_allocate_func) ((m+1) * sizeof(mpz_t));
+  ptoj = (mpz_t*) (*__gmp_allocate_func) ((m+1) * sizeof(mpz_t)); /* ptoj[i] = mantissa^(2^i) */
 #ifdef A
-  T = (mpz_t*) (*_mp_allocate_func) ((m+1) * sizeof(mpz_t));
+  T = (mpz_t*) (*__gmp_allocate_func) ((m+1) * sizeof(mpz_t));
 #endif
 #ifdef R_IS_RATIONAL
-    qtoj = (mpz_t*) (*_mp_allocate_func) ((m+1) * sizeof(mpz_t)); 
+    qtoj = (mpz_t*) (*__gmp_allocate_func) ((m+1) * sizeof(mpz_t)); 
 #endif
-  if ((P == NULL) || (S == NULL) || (ptoj == NULL)
-#ifdef A
-      || (T == NULL)
-#endif
-#ifdef R_IS_RATIONAL
-      || (qtoj == NULL)
-#endif
-      ) {
-    fprintf (stderr, "Error in mpfr_generic: no more memory available\n");
-    exit (1);
-  }
   for (i=0;i<=m;i++) { mpz_init(P[i]); mpz_init(S[i]); mpz_init(ptoj[i]);
 #ifdef R_IS_RATIONAL
   mpz_init(qtoj[i]);
@@ -191,19 +180,18 @@ GENERIC (y, p, r, m)
 	}
 
   mpz_tdiv_q(S[0], S[0], P[0]);
-  mpfr_set_z(y,S[0], GMP_RNDD);
+  mpfr_set_z(y, S[0], GMP_RNDD);
   MPFR_EXP(y) += expo; 
 
 #ifdef R_IS_RATIONAL
   /* division exacte */
   mpz_div_ui(qtoj[m], qtoj[m], r);
-  i = (MPFR_PREC(y));
-  mpfr_init2(tmp,i);
+  mpfr_init2(tmp, MPFR_PREC(y));
   mpfr_set_z(tmp, qtoj[m] , GMP_RNDD);
-  mpfr_div(y, y, tmp,GMP_RNDD);  
+  mpfr_div(y, y, tmp, GMP_RNDD);
   mpfr_clear(tmp);
 #else
-  mpfr_div_2exp(y, y, r*(i-1),GMP_RNDN); 
+  mpfr_div_2exp(y, y, r*(i-1), GMP_RNDN);
 #endif
   for (i=0;i<=m;i++) { mpz_clear(P[i]); mpz_clear(S[i]); mpz_clear(ptoj[i]); 
 #ifdef R_IS_RATIONAL
@@ -213,14 +201,14 @@ GENERIC (y, p, r, m)
   mpz_clear(T[i]);
 #endif
   }
-  (*_mp_free_func) (P, (m+1) * sizeof(mpz_t));
-  (*_mp_free_func) (S, (m+1) * sizeof(mpz_t));
-  (*_mp_free_func) (ptoj, (m+1) * sizeof(mpz_t));
+  (*__gmp_free_func) (P, (m+1) * sizeof(mpz_t));
+  (*__gmp_free_func) (S, (m+1) * sizeof(mpz_t));
+  (*__gmp_free_func) (ptoj, (m+1) * sizeof(mpz_t));
 #ifdef R_IS_RATIONAL
-  (*_mp_free_func) (qtoj, (m+1) * sizeof(mpz_t));
+  (*__gmp_free_func) (qtoj, (m+1) * sizeof(mpz_t));
 #endif
 #ifdef A
-  (*_mp_free_func) (T, (m+1) * sizeof(mpz_t));
+  (*__gmp_free_func) (T, (m+1) * sizeof(mpz_t));
 #endif
   return 0;
 }
diff --git a/mpfr/get_str.c b/mpfr/get_str.c
index 9c1391c19..4dc9a65ac 100644
--- a/mpfr/get_str.c
+++ b/mpfr/get_str.c
@@ -37,11 +37,12 @@ MA 02111-1307, USA. */
 
   For op = 3.1416 we get str = "31416" and expptr=1.
  */
+char*
 #if __STDC__
-char *mpfr_get_str(char *str, mp_exp_t *expptr, int base, size_t n,
+mpfr_get_str (char *str, mp_exp_t *expptr, int base, size_t n,
 		  mpfr_srcptr op, mp_rnd_t rnd_mode)
 #else
-char *mpfr_get_str(str, expptr, base, n, op, rnd_mode)
+mpfr_get_str (str, expptr, base, n, op, rnd_mode)
      char *str;
      mp_exp_t *expptr;
      int base;
@@ -50,8 +51,13 @@ char *mpfr_get_str(str, expptr, base, n, op, rnd_mode)
      mp_rnd_t rnd_mode;
 #endif
 {
-  double d; long e, q, div, p, err, prec, sh; mpfr_t a, b; mpz_t bz;
-  char *str0=NULL; mp_rnd_t rnd1; int f, pow2, ok=0, neg;
+  double d;
+  long e, q, div, p, err, prec, sh;
+  mpfr_t a, b;
+  mpz_t bz;
+  char *str0=NULL;
+  mp_rnd_t rnd1;
+  int f, pow2, ok=0, neg, str_is_null=0;
 
   if (base<2 || 36<base) {
     fprintf(stderr, "Error: too small or too large base in mpfr_get_str: %d\n",
@@ -62,33 +68,23 @@ char *mpfr_get_str(str, expptr, base, n, op, rnd_mode)
   neg = (MPFR_SIGN(op)<0) ? 1 : 0;
 
   if (MPFR_IS_INF(op)) { 
-    if (str == NULL) {
-      str = (*_mp_allocate_func)(neg + 4);
-      if (str == NULL) {
-	fprintf (stderr, "Error in mpfr_get_str: no more memory available\n");
-	exit (1);
-      }
-    }
-    str0 = str; 
+    if (str == NULL)
+      str = (*__gmp_allocate_func)(neg + 4);
+    str0 = str;
     if (neg) { *str++ = '-'; }
-    *str++ = 'I'; *str++ = 'n'; *str++ = 'f'; *str='\0'; 
-    return str0; 
+    *str++ = 'I'; *str++ = 'n'; *str++ = 'f'; *str='\0';
+    return str0; /* strlen(str0) = neg + 3 */
   }
 
   if (!MPFR_NOTZERO(op)) {
-    if (str == NULL) {
-      str = (*_mp_allocate_func)(neg + n + 1);
-      if (str == NULL) {
-	fprintf (stderr, "Error in mpfr_get_str: no more memory available\n");
-	exit (1);
-      }
-    }
+    if (str == NULL)
+      str = (*__gmp_allocate_func)(neg + n + 1);
     str0 = str;
-    if (MPFR_SIGN(op)<0) *str++ = '-';
+    if (neg) *str++ = '-';
     for (f=0;f<n;f++) *str++ = '0';
     *str++ = '\0';
     *expptr = 1;
-    return str0;
+    return str0; /* strlen(str0) = neg + n */
   }
 
   count_leading_zeros(pow2, (mp_limb_t)base); 
@@ -205,11 +201,7 @@ char *mpfr_get_str(str, expptr, base, n, op, rnd_mode)
   /* computes the number of characters needed */
   q = neg + n + 2; /* n+1 may not be enough for 100000... */
   if (str == NULL) {
-    str0 = str = (*_mp_allocate_func)(q);
-    if (str == NULL) {
-      fprintf (stderr, "Error in mpfr_get_str: no more memory available\n");
-      exit (1);
-    }
+    str0 = str = (*__gmp_allocate_func) (q);
   }
   if (neg) *str++='-';
   mpz_get_str(str, base, bz); /* n digits of mantissa */
@@ -224,6 +216,12 @@ char *mpfr_get_str(str, expptr, base, n, op, rnd_mode)
   }
   *expptr = f;
   mpfr_clear(a); mpfr_clear(b); mpz_clear(bz);
+
+  /* if the given string was null, ensure we return a block which is exactly
+     strlen(str)+1 bytes long (useful for __gmp_free_func and the C++ wrapper)
+  */
+  if (str_is_null && ((strlen(str0) + 1) != q))
+    str0 = (char*) (*__gmp_reallocate_func) (str0, q, strlen(str0) + 1);
+
   return str0;
 }
-
diff --git a/mpfr/init.c b/mpfr/init.c
index 7f95c4159..d2474aa78 100644
--- a/mpfr/init.c
+++ b/mpfr/init.c
@@ -43,11 +43,7 @@ mpfr_init2 (x, p)
   xsize = (p - 1)/BITS_PER_MP_LIMB + 1; 
 
   MPFR_PREC(x) = p;
-  MPFR_MANT(x) = (mp_ptr) (*_mp_allocate_func) (xsize * BYTES_PER_MP_LIMB);
-  if (MPFR_MANT(x) == NULL) {
-    fprintf (stderr, "Error in mpfr_init2: no more memory available\n");
-    exit (1);
-  }
+  MPFR_MANT(x) = (mp_ptr) (*__gmp_allocate_func) (xsize * BYTES_PER_MP_LIMB);
   MPFR_SIZE(x) = xsize;
   MPFR_CLEAR_FLAGS(x); 
   MPFR_SET_ZERO(x); /* initializes to zero */
@@ -62,5 +58,5 @@ mpfr_init (x)
      mpfr_ptr x;
 #endif
 {
-  mpfr_init2(x, __gmp_default_fp_bit_precision);
+  mpfr_init2(x, __mpfr_default_fp_bit_precision);
 }
diff --git a/mpfr/inp_str.c b/mpfr/inp_str.c
index 3a6d847f9..422e5f797 100644
--- a/mpfr/inp_str.c
+++ b/mpfr/inp_str.c
@@ -49,11 +49,7 @@ mpfr_inp_str (rop, stream, base, rnd_mode)
     stream = stdin;
 
   alloc_size = 100;
-  str = (char *) (*_mp_allocate_func) (alloc_size);
-  if (str == NULL) {
-    fprintf (stderr, "Error in mpfr_inp_str: no more memory available\n");
-    exit (1);
-  }
+  str = (char *) (*__gmp_allocate_func) (alloc_size);
   str_size = 0;
   nread = 0;
 
@@ -71,11 +67,7 @@ mpfr_inp_str (rop, stream, base, rnd_mode)
 	{
 	  size_t old_alloc_size = alloc_size;
 	  alloc_size = alloc_size * 3 / 2;
-	  str = (char *) (*_mp_reallocate_func) (str, old_alloc_size, alloc_size);
-	  if (str == NULL) {
-	    fprintf (stderr, "Error in mpfr_inp_str: no more memory available\n");
-	    exit (1);
-	  }
+	  str = (char *) (*__gmp_reallocate_func) (str, old_alloc_size, alloc_size);
 	}
       if (c == EOF || isspace (c))
 	break;
@@ -88,11 +80,7 @@ mpfr_inp_str (rop, stream, base, rnd_mode)
     {
       size_t old_alloc_size = alloc_size;
       alloc_size = alloc_size * 3 / 2;
-      str = (char *) (*_mp_reallocate_func) (str, old_alloc_size, alloc_size);
-      if (str == NULL) {
-	fprintf (stderr, "Error in mpfr_inp_str: no more memory available\n");
-	exit (1);
-      }
+      str = (char *) (*__gmp_reallocate_func) (str, old_alloc_size, alloc_size);
     }
   str[str_size] = 0;
 
@@ -100,6 +88,6 @@ mpfr_inp_str (rop, stream, base, rnd_mode)
   if (retval == -1)
     return 0;			/* error */
 
-  (*_mp_free_func) (str, alloc_size);
+  (*__gmp_free_func) (str, alloc_size);
   return str_size + nread;
 }
diff --git a/mpfr/log.c b/mpfr/log.c
index 523d19870..726a7bd5b 100644
--- a/mpfr/log.c
+++ b/mpfr/log.c
@@ -148,8 +148,8 @@ mpfr_log (r, a, rnd_mode)
     /* we have 7 ulps of error from the above roundings,
        4 ulps from the 4/s^2 second order term,
        plus the cancelled bits */
-    if (mpfr_can_round(cst,p-cancel-4,GMP_RNDN,rnd_mode,q)==1) {
-      mpfr_set(r,cst,rnd_mode);
+    if (mpfr_can_round (cst, p - cancel - 4, GMP_RNDN, rnd_mode, q) == 1) {
+      mpfr_set (r, cst, rnd_mode);
 #ifdef DEBUG
       printf("result="); mpfr_print_raw(r); putchar('\n');
 #endif
@@ -157,8 +157,7 @@ mpfr_log (r, a, rnd_mode)
     }
     /* else we increase the precision */
     else {
-      p += BITS_PER_MP_LIMB+cancel;
-      TMP_FREE(marker);
+      p += BITS_PER_MP_LIMB + cancel;
     }
 
     /* We clean */
diff --git a/mpfr/mpf2mpfr.h b/mpfr/mpf2mpfr.h
index 0e4704b8c..7b7336a22 100644
--- a/mpfr/mpf2mpfr.h
+++ b/mpfr/mpf2mpfr.h
@@ -1,86 +1,6 @@
 /* types */
 #define mpf_t mpfr_t
 
-
-#undef mpf_ceil
-#undef mpf_clear
-#undef mpf_cmp
-#undef mpf_cmp_si
-#undef mpf_cmp_ui
-#undef mpf_eq
-#undef mpf_floor
-#undef mpf_get_d
-#undef mpf_get_prec
-#undef mpf_init
-#undef mpf_init2
-#undef mpf_random2
-#undef mpf_set_default_prec
-#undef mpf_set_prec
-#undef mpf_set_prec_raw
-#undef mpf_trunc
-#undef mpf_sgn
-#undef mpf_swap
-
-/* functions which take as argument the rounding mode */
-#undef mpf_abs
-#undef mpf_add
-#undef mpf_add_ui
-#undef mpf_div
-#undef mpf_div_ui
-#undef mpf_div_2exp
-#undef mpf_dump
-#undef mpf_get_str
-#undef mpf_inp_str
-#undef mpf_set_str
-#undef mpf_init_set
-#undef mpf_init_set_d
-#undef mpf_init_set_si
-#undef mpf_init_set_str
-#undef mpf_init_set_ui
-#undef mpf_mul
-#undef mpf_mul_2exp
-#undef mpf_mul_ui
-#undef mpf_neg
-#undef mpf_out_str
-#undef mpf_pow_ui
-#undef mpf_reldiff
-#undef mpf_set
-#undef mpf_set_d
-#undef mpf_set_q
-#undef mpf_set_si
-#undef mpf_set_ui
-#undef mpf_set_z
-#undef mpf_sqrt
-#undef mpf_sqrt_ui
-#undef mpf_sub
-#undef mpf_sub_ui
-#undef mpf_ui_div
-#undef mpf_ui_sub
-#undef mpf_urandomb
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
 /* functions which don't take as argument the rounding mode */
 #define mpf_ceil mpfr_ceil
 #define mpf_clear mpfr_clear
diff --git a/mpfr/mpfr-impl.h b/mpfr/mpfr-impl.h
index 5303b8068..3e904d999 100644
--- a/mpfr/mpfr-impl.h
+++ b/mpfr/mpfr-impl.h
@@ -19,19 +19,6 @@ 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. */
 
-
-/* Fudges */
-#define _mp_allocate_func      __gmp_allocate_func
-#define _mp_reallocate_func    __gmp_reallocate_func
-#define _mp_free_func          __gmp_free_func
-#define _mp_default_allocate   __gmp_default_allocate
-#define _mp_default_reallocate __gmp_default_reallocate
-#define _mp_default_free       __gmp_default_free
-
-
-
-
-
 /* Definition of constants */
 
 #define LOG2 0.69314718055994528622 /* log(2) rounded to zero on 53 bits */
@@ -72,7 +59,6 @@ typedef union ieee_double_extract Ieee_double_extract;
 #define MPFR_NOTZERO(x) (MPFR_MANT(x)[(MPFR_PREC(x)-1)/BITS_PER_MP_LIMB])
 #define MPFR_IS_ZERO(x) ((MPFR_NOTZERO(x))==0)
 #define MPFR_SET_ZERO(x) (MPFR_MANT(x)[(MPFR_PREC(x)-1)/BITS_PER_MP_LIMB] = 0)
-#define mpfr_sgn(x) ((MPFR_NOTZERO(x)) ? MPFR_SIGN(x) : 0)
 
 /* Memory gestion */
 
diff --git a/mpfr/mpfr-test.h b/mpfr/mpfr-test.h
index ef9fbfb04..8db884ca9 100644
--- a/mpfr/mpfr-test.h
+++ b/mpfr/mpfr-test.h
@@ -29,9 +29,10 @@ MA 02111-1307, USA. */
 #define random mrand48
 #endif
 
-double drand _PROTO((void)); 
-int ulp _PROTO((double, double)); 
-double dbl _PROTO((double, int)); 
+double drand _PROTO ((void)); 
+int ulp _PROTO ((double, double)); 
+double dbl _PROTO ((double, int)); 
+double Ulp _PROTO ((double));
 
 #define MINNORM 2.2250738585072013831e-308 /* 2^(-1022), smallest normalized */
 #define MAXNORM 1.7976931348623157081e308 /* 2^(1023)*(2-2^(-52)) */
@@ -57,20 +58,34 @@ double drand ()
   return d;
 }
 
+/* returns ulp(x) for x a 'normal' double-precision number */
+double Ulp (double x)
+{
+   double y, eps;
+
+   if (x < 0) x = -x;
+
+   y = x * 2.220446049250313080847263336181640625e-16 ; /* x / 2^52 */
+
+   /* as ulp(x) <= y = x/2^52 < 2*ulp(x),
+   we have x + ulp(x) <= x + y <= x + 2*ulp(x),
+   therefore o(x + y) = x + ulp(x) or x + 2*ulp(x) */
+
+   eps =  x + y;
+   eps = eps - x; /* ulp(x) or 2*ulp(x) */
+
+   return (eps > y) ? 0.5 * eps : eps;
+}
+
 /* returns the number of ulp's between a and b */
 int ulp (double a, double b)
 {
-  double eps=1.1102230246251565404e-16; /* 2^(-53) */
   if (a==0.0) {
     if (b==0.0) return 0;
     else if (b<0.0) return 2147483647;
     else return -2147483647;
   }
-  b = (a-b)/a;
-  if (b>0)
-    return (int) floor(b/eps);
-  else
-    return (int) ceil(b/eps);
+  return (a-b)/Ulp(a);
 }
 
 /* return double m*2^e */
diff --git a/mpfr/mpfr.h b/mpfr/mpfr.h
index 948b1103e..12200338a 100644
--- a/mpfr/mpfr.h
+++ b/mpfr/mpfr.h
@@ -19,6 +19,8 @@ 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. */
 
+#ifndef __MPFR_H
+#define __MPFR_H
 #include <stdio.h>
 
 /* Definition of rounding modes */
@@ -98,12 +100,11 @@ void mpfr_set_si _PROTO ((mpfr_ptr, long, mp_rnd_t));
 void mpfr_set_ui _PROTO ((mpfr_ptr, unsigned long, mp_rnd_t)); 
 void mpfr_print_raw _PROTO ((mpfr_srcptr)); 
 void mpfr_random _PROTO ((mpfr_ptr));
-void mpfr_srandom _PROTO ((unsigned long int));
 void mpfr_random2 _PROTO ((mpfr_ptr, mp_size_t, mp_exp_t)); 
 void mpfr_urandomb _PROTO ((mpfr_ptr, gmp_randstate_t)); 
 void mpfr_clear _PROTO ((mpfr_ptr)); 
 void mpfr_set_str_raw _PROTO ((mpfr_ptr, char *));
-int mpfr_set_str _PROTO ((mpfr_ptr, char *, int, mp_rnd_t));
+int mpfr_set_str _PROTO ((mpfr_ptr, __gmp_const char *, int, mp_rnd_t));
 int mpfr_init_set_str _PROTO ((mpfr_ptr, char *, int, mp_rnd_t));
 size_t mpfr_inp_str _PROTO ((mpfr_ptr, FILE *, int, mp_rnd_t));
 char* mpfr_get_str _PROTO ((char *, mp_exp_t *, int, size_t, mpfr_srcptr, mp_rnd_t));
@@ -137,7 +138,8 @@ void mpfr_div_2exp _PROTO((mpfr_ptr, mpfr_srcptr, unsigned long int,mp_rnd_t));
 void mpfr_set_prec _PROTO((mpfr_ptr, mp_prec_t));
 void mpfr_set_prec_raw _PROTO((mpfr_ptr, mp_prec_t));
 void mpfr_set_default_prec _PROTO((mp_prec_t));
-extern mp_prec_t __gmp_default_fp_bit_precision;
+mp_prec_t mpfr_get_default_prec _PROTO((void));
+extern mp_prec_t __mpfr_default_fp_bit_precision;
 extern mp_rnd_t __gmp_default_rounding_mode;
 char * mpfr_print_rnd_mode _PROTO((mp_rnd_t)); 
 void mpfr_neg _PROTO((mpfr_ptr, mpfr_srcptr, mp_rnd_t)); 
@@ -155,6 +157,9 @@ void mpfr_swap _PROTO((mpfr_ptr, mpfr_ptr));
 void mpfr_dump _PROTO((mpfr_srcptr, mp_rnd_t));
 void mpfr_set4 _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t, int));
 int mpfr_cmp3 _PROTO ((mpfr_srcptr, mpfr_srcptr, long int));
+int mpfr_nan_p _PROTO((mpfr_srcptr));
+int mpfr_inf_p _PROTO((mpfr_srcptr));
+int mpfr_number_p _PROTO((mpfr_srcptr));
 
 #if defined (__cplusplus)
 }
@@ -165,6 +170,7 @@ int mpfr_cmp3 _PROTO ((mpfr_srcptr, mpfr_srcptr, long int));
 #define mpfr_set(a,b,r) mpfr_set4(a,b,r,MPFR_SIGN(b))
 #define mpfr_abs(a,b,r) mpfr_set4(a,b,r,1)
 #define mpfr_cmp(b, c) mpfr_cmp3(b, c, 1)
+#define mpfr_sgn(x) ((MPFR_NOTZERO(x)) ? MPFR_SIGN(x) : 0)
 
 #define mpfr_init_set_si(x, i, rnd) \
  do { mpfr_init(x); mpfr_set_si((x), (i), (rnd)); } while (0)
@@ -176,3 +182,5 @@ int mpfr_cmp3 _PROTO ((mpfr_srcptr, mpfr_srcptr, long int));
  do { mpfr_init(x); mpfr_set((x), (y), (rnd)); } while (0)
 #define mpfr_init_set_f(x, y, rnd) \
  do { mpfr_init(x); mpfr_set_f((x), (y), (rnd)); } while (0)
+
+#endif
diff --git a/mpfr/mpz_set_fr.c b/mpfr/mpz_set_fr.c
index 92f14f17d..007bf0e9e 100644
--- a/mpfr/mpz_set_fr.c
+++ b/mpfr/mpz_set_fr.c
@@ -43,15 +43,7 @@ mpz_set_fr (z, f)
   fn = 1 + (MPFR_PREC(f)-1)/BITS_PER_MP_LIMB;
 
   /* check whether allocated space for z is enough */
-  if (ALLOC(z) < fn) {
-    PTR(z) = (mp_ptr) (*_mp_reallocate_func) (PTR(z), 
-	   ABS(SIZ(z))*BYTES_PER_MP_LIMB, fn*BYTES_PER_MP_LIMB);
-    if (PTR(z) == NULL) {
-      fprintf (stderr, "Error in mpz_set_fr: no more memory available\n");
-      exit (1);
-    }
-    ALLOC(z) = fn;
-  }
+  if (ALLOC(z) < fn) MPZ_REALLOC(z,fn);
 
   sh = fn*BITS_PER_MP_LIMB - MPFR_PREC(f);
   if (sh) mpn_rshift(PTR(z), MPFR_MANT(f), fn, sh);
diff --git a/mpfr/out_str.c b/mpfr/out_str.c
index 9b7d7408b..3f7ebbbcb 100644
--- a/mpfr/out_str.c
+++ b/mpfr/out_str.c
@@ -39,7 +39,9 @@ mpfr_out_str (stream, base, n_digits, op, rnd_mode)
      mp_rnd_t rnd_mode;
 #endif
 {
-  char *s,*s0; size_t l; mp_exp_t e;
+  char *s, *s0;
+  size_t l;
+  mp_exp_t e;
 
   if (MPFR_IS_NAN(op)) { fprintf(stream, "NaN"); return 3; }
   if (!MPFR_NOTZERO(op)) { fprintf(stream, "0"); return 1; }
@@ -49,7 +51,7 @@ mpfr_out_str (stream, base, n_digits, op, rnd_mode)
       else { fprintf(stream, "-Inf"); return 4; }
     }
 
-  s = mpfr_get_str(NULL, &e, base, n_digits, op, rnd_mode);
+  s = mpfr_get_str (NULL, &e, base, n_digits, op, rnd_mode);
   /*	  TODO: maintenir le code pour les infinis dans get_str ?  */
   s0 = s;
   /* for op=3.1416 we have s = "31416" and e = 1 */
@@ -61,12 +63,9 @@ mpfr_out_str (stream, base, n_digits, op, rnd_mode)
   fputc('.', stream);       /* decimal point */
   fputs(s, stream);         /* rest of mantissa */
   if (e) {
-    fputc((base>10) ? '@' : 'e', stream); l++;
-    sprintf(s, "%ld", e);
-    l += strlen(s);
-    fprintf(stream, "%s", s);
+    l += fprintf (stream, (base <= 10 ? "e%ld" : "@%ld"), e);
   }
 
-  (*_mp_free_func)(s0, l); 
+  (*__gmp_free_func)(s0, l); 
   return l;
 }
diff --git a/mpfr/print_raw.c b/mpfr/print_raw.c
index 8d2e44212..3b75e8eb3 100644
--- a/mpfr/print_raw.c
+++ b/mpfr/print_raw.c
@@ -94,15 +94,11 @@ mpfr_print_raw(x)
 	= 17 + MPFR_ABSSIZE(x) * BITS_PER_MP_LIMB
       */
     alloc_size = 17 + MPFR_ABSSIZE(x) * BITS_PER_MP_LIMB;
-     str = (char *) (*_mp_allocate_func) (alloc_size * sizeof(char));
-     if (str == NULL) {
-       fprintf (stderr, "Error in mpfr_print_raw: no more memory available\n");
-       exit (1);
-     }
+     str = (char *) (*__gmp_allocate_func) (alloc_size * sizeof(char));
      mpfr_get_str_raw(str, x);
 
      printf("%s", str); 
-     (*_mp_free_func) (str, alloc_size * sizeof(char));
+     (*__gmp_free_func) (str, alloc_size * sizeof(char));
   }
 }
 
diff --git a/mpfr/random.c b/mpfr/random.c
index e2a50770b..64532b629 100644
--- a/mpfr/random.c
+++ b/mpfr/random.c
@@ -55,16 +55,3 @@ mpfr_random (x)
   /* cnt is the number of non significant bits in the low limb */
   xp[0] &= ~((((mp_limb_t) 1) << cnt) - 1);
 }
-
-#if HAVE_SRANDOM
-void
-#if __STDC__
-mpfr_srandom (unsigned long int seed)
-#else
-mpfr_srandom (seed)
-     unsigned long int seed;
-#endif
-{
-  srandom (seed); 
-}
-#endif
diff --git a/mpfr/round.c b/mpfr/round.c
index 06de436b9..a6de9bd50 100644
--- a/mpfr/round.c
+++ b/mpfr/round.c
@@ -113,10 +113,10 @@ mpfr_round_raw2(xp, xn, neg, rnd, prec)
 */
 int
 #if __STDC__
-mpfr_round_raw(mp_limb_t *y, mp_limb_t *xp, mp_prec_t xprec, int negative,
+mpfr_round_raw (mp_limb_t *y, mp_limb_t *xp, mp_prec_t xprec, int negative,
 	       mp_prec_t yprec, mp_rnd_t rnd_mode)
 #else
-mpfr_round_raw(y, xp, xprec, negative, yprec, rnd_mode)
+mpfr_round_raw (y, xp, xprec, negative, yprec, rnd_mode)
      mp_limb_t *y; 
      mp_limb_t *xp; 
      mp_prec_t xprec; 
@@ -142,7 +142,8 @@ mpfr_round_raw(y, xp, xprec, negative, yprec, rnd_mode)
   /* precision is larger than the size of x. No rounding is necessary. 
      Just add zeroes at the end */
   if (xsize < nw) { 
-    MPN_COPY(y + nw - xsize, xp, xsize);
+    /* if y=xp, maybe an overlap: MPN_COPY_DECR is ok when src <= dst */
+    MPN_COPY_DECR(y + nw - xsize, xp, xsize);
     MPN_ZERO(y, nw - xsize); /* PZ 27 May 99 */
     y[0] &= mask;
     return 0; 
@@ -156,8 +157,8 @@ mpfr_round_raw(y, xp, xprec, negative, yprec, rnd_mode)
       else
 	carry = mpn_add_1(y, xp + xsize - nw, nw, 1); 
     }
-
-  else MPN_COPY(y, xp + xsize - nw, nw);
+  else /* now xsize >= nw */
+    MPN_COPY_INCR(y, xp + xsize - nw, nw);
 
   y[0] &= mask;
   return carry; 
@@ -165,9 +166,9 @@ mpfr_round_raw(y, xp, xprec, negative, yprec, rnd_mode)
 
 void
 #if __STDC__
-mpfr_round(mpfr_ptr x, mp_rnd_t rnd_mode, mp_prec_t prec)
+mpfr_round (mpfr_ptr x, mp_rnd_t rnd_mode, mp_prec_t prec)
 #else
-mpfr_round(x, rnd_mode, prec)
+mpfr_round (x, rnd_mode, prec)
      mpfr_ptr x; 
      mp_rnd_t rnd_mode; 
      mp_prec_t prec; 
@@ -183,12 +184,8 @@ mpfr_round(x, rnd_mode, prec)
 
   /* check if x has enough allocated space for the mantissa */
   if (nw > MPFR_ABSSIZE(x)) {
-    MPFR_MANT(x) = (mp_ptr) (*_mp_reallocate_func) 
+    MPFR_MANT(x) = (mp_ptr) (*__gmp_reallocate_func) 
       (MPFR_MANT(x), MPFR_ABSSIZE(x)*BYTES_PER_MP_LIMB, nw * BYTES_PER_MP_LIMB);
-    if (MPFR_MANT(x) == NULL) {
-      fprintf (stderr, "Error in mpfr_round: no more memory available\n");
-      exit (1);
-    }
     MPFR_SIZE(x) = nw; /* new number of allocated limbs */
   }
 
diff --git a/mpfr/set_d.c b/mpfr/set_d.c
index 3084e71cb..113066ce4 100644
--- a/mpfr/set_d.c
+++ b/mpfr/set_d.c
@@ -338,9 +338,19 @@ mpfr_get_d2(src, e)
   /* Accumulate the limbs from less significant to most significant
      otherwise due to rounding we may accumulate several ulps,
      especially in rounding towards -/+infinity. */
-  for (i = n_limbs_to_use; i>=1; i--)
+  for (i = n_limbs_to_use; i>=1; i--) {
+#if (BITS_PER_MP_LIMB == 32)
     res = res / MP_BASE_AS_DOUBLE +
       ((negative) ? -(double)qp[size - i] : qp[size - i]);
+#else
+#if (BITS_PER_MP_LIMB == 64)
+    q = qp[size - i] & (mp_limb_t) 4294967295;
+    res = res / MP_BASE_AS_DOUBLE + ((negative) ? -(double)q : q);
+    q = qp[size - i] - q;
+    res = res + ((negative) ? -(double)q : q);
+#endif /* BITS_PER_MP_LIMB == 64 */
+#endif /* BITS_PER_MP_LIMB == 32 */
+  }
   res = __mpfr_scale2 (res, e - BITS_PER_MP_LIMB); 
 
   return res;
diff --git a/mpfr/set_dfl_prec.c b/mpfr/set_dfl_prec.c
index f1ff78a08..ecfd11ac9 100644
--- a/mpfr/set_dfl_prec.c
+++ b/mpfr/set_dfl_prec.c
@@ -1,6 +1,6 @@
-/* mpfr_set_default_prec -- set the default precision
+/* mpfr_set_default_prec, mpfr_get_default_prec -- set/get default precision
 
-Copyright (C) 1999 Free Software Foundation.
+Copyright (C) 1999, 2000 Free Software Foundation.
 
 This file is part of the MPFR Library.
 
@@ -24,7 +24,7 @@ MA 02111-1307, USA. */
 #include "gmp-impl.h"
 
 /* default is 53 bits */
-mp_prec_t __gmp_default_fp_bit_precision = 53;
+mp_prec_t __mpfr_default_fp_bit_precision = 53;
 
 void
 #if __STDC__
@@ -34,5 +34,11 @@ mpfr_set_default_prec (prec_in_bits)
      mp_prec_t prec_in_bits;
 #endif
 {
-  __gmp_default_fp_bit_precision = prec_in_bits;
+  __mpfr_default_fp_bit_precision = prec_in_bits;
+}
+
+mp_prec_t
+mpfr_get_default_prec ()
+{
+  return __mpfr_default_fp_bit_precision;
 }
diff --git a/mpfr/set_prec.c b/mpfr/set_prec.c
index 17adb5bb2..3f52e5b17 100644
--- a/mpfr/set_prec.c
+++ b/mpfr/set_prec.c
@@ -43,12 +43,8 @@ mpfr_set_prec (x, p)
   xsize = (p - 1)/BITS_PER_MP_LIMB + 1; /* new limb size */
 
   if (xsize > MPFR_ABSSIZE(x)) {
-    MPFR_MANT(x) = (mp_ptr) (*_mp_reallocate_func) (MPFR_MANT(x), 
+    MPFR_MANT(x) = (mp_ptr) (*__gmp_reallocate_func) (MPFR_MANT(x), 
 		 MPFR_ABSSIZE(x)*BYTES_PER_MP_LIMB, xsize * BYTES_PER_MP_LIMB);
-    if (MPFR_MANT(x) == NULL) {
-      fprintf (stderr, "Error in mpfr_set_prec: no more memory available\n");
-      exit (1);
-    }
     MPFR_SIZE(x) = xsize; /* new number of allocated limbs */
   }
 
diff --git a/mpfr/set_si.c b/mpfr/set_si.c
index 9594682fe..20abb0e7f 100644
--- a/mpfr/set_si.c
+++ b/mpfr/set_si.c
@@ -1,6 +1,6 @@
-/* mpfr_set_si, mpfr_set_ui -- set a MPFR number from a machine integer
+/* mpfr_set_si -- set a MPFR number from a machine signed integer
 
-Copyright (C) 1999 Free Software Foundation.
+Copyright (C) 1999-2001 Free Software Foundation.
 
 This file is part of the MPFR Library.
 
@@ -41,7 +41,7 @@ mpfr_set_si (x, i, rnd_mode)
   MPFR_CLEAR_FLAGS(x);
   if (i==0) { MPFR_SET_ZERO(x); return; }
   xn = (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB;
-  ai = ABS(i); 
+  ai = (unsigned long) ABS(i); 
 
   count_leading_zeros(cnt, ai); 
 
@@ -67,44 +67,3 @@ mpfr_set_si (x, i, rnd_mode)
 
   return; 
 }
-
-void
-#if __STDC__
-mpfr_set_ui (mpfr_ptr x, unsigned long int i, mp_rnd_t rnd_mode)
-#else
-mpfr_set_ui (x, i, rnd_mode)
-     mpfr_ptr x;
-     unsigned long int i;
-     mp_rnd_t rnd_mode;
-#endif  
-{
-  unsigned int xn, cnt; mp_limb_t *xp;
-
-  MPFR_CLEAR_FLAGS(x);
-  if (i==0) { MPFR_SET_ZERO(x); return; }
-  xn = (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB;
-  count_leading_zeros(cnt, (mp_limb_t) i); 
-
-  xp = MPFR_MANT(x);
-  xp[xn] = ((mp_limb_t) i) << cnt; 
-  /* don't forget to put zero in lower limbs */
-  MPN_ZERO(xp, xn);
-
-  MPFR_EXP(x) = BITS_PER_MP_LIMB - cnt;
-
-  /* round if MPFR_PREC(x) smaller than length of i */
-  if (MPFR_PREC(x) < BITS_PER_MP_LIMB-cnt) {
-    cnt = mpfr_round_raw(xp+xn, xp+xn, BITS_PER_MP_LIMB-cnt, 0, MPFR_PREC(x), 
-			 rnd_mode);
-    if (cnt) { /* special case 1.000...000 */
-      MPFR_EXP(x)++;
-      xp[xn] = ((mp_limb_t) 1) << (BITS_PER_MP_LIMB-1);
-    }
-  }
-
-  /* warning: don't change the precision of x! */
-  if (MPFR_SIGN(x) < 0) MPFR_CHANGE_SIGN(x);
-
-  return; 
-}
-
diff --git a/mpfr/set_str.c b/mpfr/set_str.c
index 0c0b029cf..b97081765 100644
--- a/mpfr/set_str.c
+++ b/mpfr/set_str.c
@@ -35,11 +35,11 @@ MA 02111-1307, USA. */
 
 int
 #if __STDC__
-mpfr_set_str(mpfr_ptr x, char *str, int base, mp_rnd_t rnd_mode)
+mpfr_set_str (mpfr_ptr x, __gmp_const char *str, int base, mp_rnd_t rnd_mode)
 #else
-mpfr_set_str(x, str, base, rnd_mode)
+mpfr_set_str (x, str, base, rnd_mode)
      mpfr_ptr x;
-     char *str;
+     __gmp_const char *str;
      int base;
      mp_rnd_t rnd_mode;
 #endif
@@ -98,7 +98,7 @@ mpfr_set_str(x, str, base, rnd_mode)
   }
   else {
     expn = -k;
-    endptr = str;
+    endptr = (char*) str;
   }
 
   /* the number is mantissa*base^expn */
diff --git a/mpfr/set_str_raw.c b/mpfr/set_str_raw.c
index e00645cef..4561daa46 100644
--- a/mpfr/set_str_raw.c
+++ b/mpfr/set_str_raw.c
@@ -1,6 +1,6 @@
 /* mpfr_set_str_raw -- set a floating-point number from a binary string
 
-Copyright (C) 1999 Free Software Foundation.
+Copyright (C) 1999-2001 Free Software Foundation.
 
 This file is part of the MPFR Library.
 
@@ -52,12 +52,7 @@ mpfr_set_str_raw (x, str)
   xp = MPFR_MANT(x);
   xsize = 1 + (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB;
   alloc = (strlen(str)+1) * sizeof(char);
-  str0 = str2 = (char *) (*_mp_allocate_func) (alloc);
-
-  if (str0 == NULL) {
-    fprintf (stderr, "Error in mpfr_set_str_raw: no more memory available\n");
-    exit (1);
-  }
+  str0 = str2 = (char *) (*__gmp_allocate_func) (alloc);
 
   if (*str == '-') { negative = 1; str++; }
   else if (*str == '+') str++;
@@ -103,35 +98,42 @@ mpfr_set_str_raw (x, str)
   else expn=k;
 
   endstr2 = str2;
-  *str2 = (char) 0; /* end of string */
-  l = (strlen(str0) - 1) / BITS_PER_MP_LIMB + 1; str2 = str0;
-  if (l > xsize) {
-    fprintf (stderr, "Error: mantissa larger than precision of destination variable in mpfr_set_str_raw\n");
-    exit (1);
+  l = endstr2 - str0; /* length of mantissa */
+  if (l == 0) { /* input is zero */
+    MPFR_SET_ZERO(x);
   }
+  else {
+    l = (l-1) / BITS_PER_MP_LIMB + 1;
+    str2 = str0;
 
-  /* str2[0]..endstr2[-1] contains the mantissa */
-  for (k = 1; k <= l; k++)
-    {
-      j = 0; 
-      xp[xsize - k] = 0; 
-      while (str2<endstr2 && j < BITS_PER_MP_LIMB)
-	{
-	  xp[xsize - k] = (xp[xsize - k] << 1) + (*str2 - '0'); 
-	  str2++; j++; 
-	}
-      xp[xsize - k] <<= (BITS_PER_MP_LIMB - j); 
+    if (l > xsize) {
+      fprintf (stderr, "Error: mantissa larger than precision of destination variable in mpfr_set_str_raw\n");
+      exit (1);
     }
 
-  for (; k <= xsize; k++) { xp[xsize - k] = 0; }
-
-  count_leading_zeros(cnt, xp[xsize - 1]); 
-  if (cnt) mpn_lshift(xp, xp, xsize, cnt); 
-
-  MPFR_EXP(x) = expn - cnt; 
-  MPFR_SIZE(x) = xsize; if (negative) MPFR_CHANGE_SIGN(x);
+    /* str2[0]..endstr2[-1] contains the mantissa */
+    for (k = 1; k <= l; k++)
+      {
+	j = 0; 
+	xp[xsize - k] = 0; 
+	while (str2<endstr2 && j < BITS_PER_MP_LIMB)
+	  {
+	    xp[xsize - k] = (xp[xsize - k] << 1) + (*str2 - '0'); 
+	    str2++; j++; 
+	  }
+	xp[xsize - k] <<= (BITS_PER_MP_LIMB - j); 
+      }
+
+    for (; k <= xsize; k++) { xp[xsize - k] = 0; }
+
+    count_leading_zeros(cnt, xp[xsize - 1]); 
+    if (cnt) mpn_lshift(xp, xp, xsize, cnt); 
+
+    MPFR_EXP(x) = expn - cnt; 
+    MPFR_SIZE(x) = xsize; if (negative) MPFR_CHANGE_SIGN(x);
+  }
 
-  (*_mp_free_func) (str0, alloc);
+  (*__gmp_free_func) (str0, alloc);
   
   /* May change to take into account : syntax errors, overflow in exponent, 
      string truncated because of size of x */
diff --git a/mpfr/sin_cos.c b/mpfr/sin_cos.c
index fc71eb2f0..ff1300f09 100644
--- a/mpfr/sin_cos.c
+++ b/mpfr/sin_cos.c
@@ -71,7 +71,7 @@ mpfr_sin_cos (sinus, cosinus, x, rnd_mode)
   mpfr_t tmp;
   mpfr_t inter;
   int ttt;
-  int twopoweri;
+  int twopoweri, expo;
   int Prec;
   int loop;
   int prec_x;
@@ -83,6 +83,7 @@ mpfr_sin_cos (sinus, cosinus, x, rnd_mode)
   int logn;
   int tmp_factor;
   int tmpi;
+  TMP_DECL (marker);
 
   if (sinus == cosinus) {
     fprintf (stderr, "Error in mpfr_sin_cos: 1st and 2nd operands must be different\n");
@@ -90,22 +91,45 @@ mpfr_sin_cos (sinus, cosinus, x, rnd_mode)
   }
 
   if (MPFR_IS_NAN(x) || MPFR_IS_INF(x)) {
-    MPFR_SET_NAN(sinus);
-    MPFR_SET_NAN(cosinus);
+    if (sinus != NULL) MPFR_SET_NAN(sinus);
+    if (cosinus != NULL) MPFR_SET_NAN(cosinus);
     return 1; /* inexact */
   }
 
   if (!MPFR_NOTZERO(x)) { 
-    mpfr_set_ui(sinus, 0, GMP_RNDN); 
-    mpfr_set_ui(cosinus, 1, GMP_RNDN); 
+    if (sinus != NULL) mpfr_set_ui (sinus, 0, GMP_RNDN); 
+    if (cosinus != NULL) mpfr_set_ui (cosinus, 1, GMP_RNDN); 
     return 0; /* exact results */
   }
 
+  TMP_MARK (marker);
+  /* allow sinus or cosinus to be NULL */
+  if (sinus == NULL) { mp_size_t s;
+    if (cosinus == NULL) {
+      fprintf (stderr, "Error in mpfr_sin_cos: 1st and 2nd operands cannot be NULL simultaneously\n");
+      exit (1);
+    }
+    s = 1 + (MPFR_PREC(cosinus) - 1) / BITS_PER_MP_LIMB;
+    sinus = TMP_ALLOC (sizeof(mpfr_t));
+    MPFR_MANT(sinus) = (mp_ptr) TMP_ALLOC(s*BYTES_PER_MP_LIMB);
+    MPFR_PREC(sinus) = MPFR_PREC(cosinus);
+    MPFR_SIZE(sinus) = s;
+    MPFR_EXP(sinus) = 0;
+  }
+  else if (cosinus == NULL) { mp_size_t s;
+    s = 1 + (MPFR_PREC(sinus) - 1) / BITS_PER_MP_LIMB;
+    cosinus = TMP_ALLOC (sizeof(mpfr_t));
+    MPFR_MANT(cosinus) = (mp_ptr) TMP_ALLOC(s*BYTES_PER_MP_LIMB);
+    MPFR_PREC(cosinus) = MPFR_PREC(sinus);
+    MPFR_SIZE(cosinus) = s;
+    MPFR_EXP(cosinus) = 0;
+  }
+
   prec_x = _mpfr_ceil_log2 ((double) MPFR_PREC(x) / BITS_PER_MP_LIMB);
   ttt = MPFR_EXP(x);
   mpfr_init2(x_copy, MPFR_PREC(x));
   mpfr_set(x_copy,x,GMP_RNDD);
-  mpz_init(square);
+  mpz_init (square);
   /* on fait le shift pour que le nombre soit inferieur a 1 */
   if (ttt > 0) 
     {
@@ -135,25 +159,31 @@ mpfr_sin_cos (sinus, cosinus, x, rnd_mode)
     mpfr_set_ui(tmp_cos,1,GMP_RNDN);
     twopoweri = BITS_PER_MP_LIMB;   
     if (k <= prec_x) iter = k; else iter= prec_x;
-    for(i = 0; i <= iter; i++){
+    for(i = 0; i <= iter; i++) {
       mpfr_extract (uk, x_copy, i);
+      expo = twopoweri - ttt;
 	if (i)
 	  {
+	    if (twopoweri > target_prec) {
+	      /* truncate the last uk, i.e. bits twopoweri/2 to twopoweri */
+	      mpz_div_2exp (uk, uk, twopoweri - target_prec);
+	      expo -= twopoweri - target_prec;
+	    }
 	    mpz_mul (square, uk, uk);
 	    mpz_neg (square, square);
-	    mpfr_sin_aux (t_sin, square, 2*(twopoweri - ttt) + 2, k - i + 1);
-	    mpfr_cos_aux (t_cos, square, 2*(twopoweri - ttt) + 2, k - i + 1);
+	    mpfr_sin_aux (t_sin, square, 2*expo + 2, k - i + 1);
+	    mpfr_cos_aux (t_cos, square, 2*expo + 2, k - i + 1);
 	    mpfr_set_z (tmp, uk, GMP_RNDD);
 	    mpfr_mul (t_sin, t_sin, tmp, GMP_RNDD);
-	    mpfr_div_2exp (t_sin, t_sin, twopoweri - ttt, GMP_RNDD);
+	    mpfr_div_2exp (t_sin, t_sin, MIN(twopoweri, target_prec) - ttt, GMP_RNDD);
 	  }
 	else
 	  {
 	    /* cas particulier : on est oblige de faire les calculs avec x/2^. 
 	       puis elever au carre (plus rapide) */    
 	    mpz_set (square, uk);
-	    mpz_mul(square, square, square);
-	    mpz_neg(square, square);
+	    mpz_mul (square, square, square);
+	    mpz_neg (square, square);
 	    /* pour l'instant, shift = 0 ... */
 	    /* ATTENTION, IL FAUT AUSSI MULTIPLIER LE DENOMINATEUR */
 	    mpfr_sin_aux(t_sin,square, 2*(shift + twopoweri - ttt) + 2, k+1);
@@ -239,8 +269,7 @@ mpfr_sin_cos (sinus, cosinus, x, rnd_mode)
   mpz_clear (uk);
   mpz_clear (square);
   mpfr_clear (x_copy);
+  TMP_FREE (marker);
+
   return 1; /* inexact result */
 } 
-
-
-
diff --git a/mpfr/sqrtrem.c b/mpfr/sqrtrem.c
index 6e16ca598..35705e0bc 100644
--- a/mpfr/sqrtrem.c
+++ b/mpfr/sqrtrem.c
@@ -1,7 +1,7 @@
 /* mpn_sqrtrem_new -- integer square root with remainder
    (should be directly integrated in a future release of GNU MP)
 
-Copyright (C) 2000 Free Software Foundation.
+Copyright (C) 1999, 2000 Free Software Foundation.
 
 This file is part of the MPFR Library.
 
@@ -25,10 +25,10 @@ MA 02111-1307, USA. */
 #include "gmp-impl.h"
 #include "longlong.h"
 
-mp_size_t mpn_sqrtrem1 _PROTO ((mp_ptr, mp_ptr, mp_srcptr));
-mp_limb_t mpn_sqrtrem2 _PROTO ((mp_ptr, mp_ptr, mp_srcptr));
-mp_limb_t mpn_dq_sqrtrem _PROTO ((mp_ptr, mp_ptr, mp_size_t));
-mp_size_t mpn_sqrtrem_new _PROTO ((mp_ptr, mp_ptr, mp_srcptr, mp_size_t));
+mp_size_t mpn_sqrtrem1 (mp_ptr, mp_ptr, mp_srcptr);
+mp_limb_t mpn_sqrtrem2 (mp_ptr, mp_ptr, mp_srcptr);
+mp_limb_t mpn_dq_sqrtrem (mp_ptr, mp_ptr, mp_size_t);
+mp_size_t mpn_sqrtrem_new (mp_ptr, mp_ptr, mp_srcptr, mp_size_t);
 
 /* table generated by isqrt(2^8*i) $ i=64..255 in MuPAD */
 static const unsigned char approx_tab[192] = {
@@ -155,7 +155,7 @@ mpn_dq_sqrtrem (mp_ptr sp, mp_ptr np, mp_size_t n)
   sp[l-1] |= q << (BITS_PER_MP_LIMB - 1);
   q >>= 1;
   if (c) c = mpn_add_n (np + l, np + l, sp + l, h);
-  mpn_mul_n (np + n, sp, sp, l);
+  mpn_sqr_n (np + n, sp, l);
   b = q + mpn_sub_n (np, np, np + n, 2 * l);
   c -= (l == h) ? b : mpn_sub_1 (np + 2 * l, np + 2 * l, 1, b);
   q = mpn_add_1 (sp + l, sp + l, h, q);
@@ -173,21 +173,19 @@ mpn_dq_sqrtrem (mp_ptr sp, mp_ptr np, mp_size_t n)
 mp_size_t
 mpn_sqrtrem_new (mp_ptr sp, mp_ptr rp, mp_srcptr np, mp_size_t nn)
 {
-  mp_limb_t *tp, s0[1], cc;
-  int c, tn;
-  mp_size_t rn;
+  mp_limb_t *tp, s0[1], cc, high, rl;
+  int c;
+  mp_size_t rn, tn;
   TMP_DECL (marker);
 
-#ifdef DEBUG
-  if (np[nn-1] == 0) {
-    fprintf (stderr, "Error in mpn_sqrtrem: most significant limb is zero\n");
-    exit(1);
-  }
-#endif
+  /* If OP is zero, both results are zero.  */
+  if (nn == 0)
+    return 0;
 
-  count_leading_zeros(c, np[nn-1]);
-
-  if (nn == 1 && c == 0) return mpn_sqrtrem1 (sp, rp, np);
+  high = np[nn-1];
+  if (nn == 1 && (high & MP_LIMB_T_HIGHBIT))
+    return mpn_sqrtrem1 (sp, rp, np);
+  count_leading_zeros(c, high);
 
   c = c / 2; /* we have to shift left by 2c bits to normalize {np, nn} */
   tn = (nn+1) / 2; /* 2*tn is the smallest even integer >= nn */
@@ -198,21 +196,20 @@ mpn_sqrtrem_new (mp_ptr sp, mp_ptr rp, mp_srcptr np, mp_size_t nn)
     tp[0] = 0; /* needed only when 2*tn > nn, but saves a test */
     if (c) mpn_lshift(tp + 2*tn - nn, np, nn, 2 * c);
     else MPN_COPY (tp + 2*tn - nn, np, nn);
-    rn = mpn_dq_sqrtrem (sp, tp, tn);
+    rl = mpn_dq_sqrtrem (sp, tp, tn);
     /* we have 2^(2k)*N = S^2 + R where k = c + (2tn-nn)*BITS_PER_MP_LIMB/2,
        thus 2^(2k)*N = (S-s0)^2 + 2*S*s0 - s0^2 + R where s0=S mod 2^k */
     c += (nn % 2) * BITS_PER_MP_LIMB / 2; /* c now represents k */
     s0[0] = sp[0] & (((mp_limb_t) 1 << c) - 1); /* S mod 2^k */
-    rn += mpn_addmul_1 (tp, sp, tn, 2 * s0[0]); /* R = R + 2*s0*S */
+    rl += mpn_addmul_1 (tp, sp, tn, 2 * s0[0]); /* R = R + 2*s0*S */
     cc = mpn_submul_1 (tp, s0, 1, s0[0]);
-    rn -= (tn > 1) ? mpn_sub_1(tp + 1, tp + 1, tn - 1, cc) : cc;
+    rl -= (tn > 1) ? mpn_sub_1(tp + 1, tp + 1, tn - 1, cc) : cc;
     mpn_rshift (sp, sp, tn, c);
-    tp[tn] = rn;
+    tp[tn] = rl;
     if (rp == NULL) rp = tp;
     c = c << 1;
     if (c < BITS_PER_MP_LIMB) tn++; else { tp++; c -= BITS_PER_MP_LIMB; }
-    if (c) mpn_rshift (rp, tp, tn, c);
-    else MPN_COPY(rp, tp, tn);
+    if (c) mpn_rshift (rp, tp, tn, c); else MPN_COPY(rp, tp, tn);
     rn = tn;
   }
   else {
@@ -221,9 +218,10 @@ mpn_sqrtrem_new (mp_ptr sp, mp_ptr rp, mp_srcptr np, mp_size_t nn)
     if (rp != np) MPN_COPY (rp, np, nn);
     rn = tn + (rp[tn] = mpn_dq_sqrtrem (sp, rp, tn));
   }
-  while (rn && rp[rn-1]==0) rn--;
-  TMP_FREE (marker);
+  
+  MPN_NORMALIZE (rp, rn);
 
+  TMP_FREE (marker);
   return rn;
 }
 
diff --git a/mpfr/tests/tlog.c b/mpfr/tests/tlog.c
index 95b27610a..7c0a32d5f 100644
--- a/mpfr/tests/tlog.c
+++ b/mpfr/tests/tlog.c
@@ -28,18 +28,25 @@ MA 02111-1307, USA. */
 #include "mpfr-test.h"
 
 double drand_log _PROTO((void)); 
-int check1 _PROTO((double, mp_rnd_t, double, int)); 
+int check1 _PROTO((double, mp_rnd_t, double, int, int)); 
 void check3 _PROTO((double, unsigned long, mp_rnd_t)); 
 void check4 _PROTO((int)); 
 void slave _PROTO((int, int)); 
 void check_worst_cases _PROTO((void)); 
 void special _PROTO((void));
 
+#if (BITS_PER_LONGINT == 32)
+#define INT32 long int
+#else
+#define INT32 int
+#endif
+
+
 double drand_log ()
 {
-  double d; long int *i;
+  double d; INT32 *i;
 
-  i = (long int*) &d;
+  i = (INT32*) &d;
   do {
     i[0] = lrand48();
     i[1] = lrand48();
@@ -48,13 +55,14 @@ double drand_log ()
   return d;
 }
 
-#define check2(a,rnd,res) check1(a,rnd,res,1)
+#define check2(a,rnd,res) check1(a,rnd,res,1,0)
 #define check(a,r) check2(a,r,0.0)
 
-int check1 (double a, mp_rnd_t rnd_mode, double res1, int ck)
+int check1 (double a, mp_rnd_t rnd_mode, double res1, int ck, int max_ulp)
 {
   mpfr_t ta, tres;
   double res2;
+  int diff=0;
   /* ck=1 iff res1 is certified correct */
 
 #ifdef TEST
@@ -69,20 +77,21 @@ int check1 (double a, mp_rnd_t rnd_mode, double res1, int ck)
   mpfr_clear(ta); mpfr_clear(tres); 
 
   if (res1!=res2 && (!isnan(res1) || !isnan(res2))) {
+      diff = ulp(res1,res2);
       if (ck) { 
 	printf("mpfr_log failed for    a=%1.20e, rnd_mode=%s\n", a,
 	       mpfr_print_rnd_mode(rnd_mode));
 	printf("correct result is        %1.20e\n mpfr_log gives          %1.20e (%d ulp)\n",res1,res2,ulp(res1,res2));
 	exit(1);
       }
-      else {
+      else if (diff>max_ulp) {
 	printf("mpfr_log differs from libm.a for a=%1.20e, rnd_mode=%s\n", a,
 	       mpfr_print_rnd_mode(rnd_mode));
 	printf(" double calculus gives %1.20e\n mpfr_log        gives %1.20e (%d ulp)\n",res1,res2,ulp(res1,res2));
       }
   }
   if (!isnan(res1) || !isnan(res2))
-    return ulp(res1,res2);
+    return diff;
   else
     return 0;
 }
@@ -101,12 +110,14 @@ void check3 (double d, unsigned long prec, mp_rnd_t rnd)
 
 void check4 (int N)
 {
-  int i, max=-1, sum=0, cur;
+  int i, max=0, sum=0, cur;
   double d;
+  mp_rnd_t rnd;
 
   for(i=0;i<N;i++) {
     d=drand_log();
-    cur=check1(d,rand() % 4,0.0,0);
+    rnd = rand() % 4;
+    cur=check1 (d, rnd, 0.0, 0, max);
     if (cur<0)
       cur = -cur;
     if (cur > max)
@@ -114,7 +125,7 @@ void check4 (int N)
     sum+=cur;
   }
   d=(double)sum / (double)N;
-  printf("max error : %i \t mean error : %f   (in ulps)\n",max,d);
+  fprintf(stderr, "max error : %i \t mean error : %f   (in ulps)\n",max,d);
 }
 
 void slave (int N, int p)
@@ -233,6 +244,14 @@ void special ()
     fprintf (stderr, "Error in mpfr_log(3) for GMP_RNDD\n");
     exit (1);
   }
+
+  /* check large precision */
+  mpfr_set_prec (x, 3322);
+  mpfr_set_prec (y, 3322);
+  mpfr_set_d (x, 3.0, GMP_RNDN);
+  mpfr_sqrt (x, x, GMP_RNDN);
+  mpfr_log (y, x, GMP_RNDN);
+
   mpfr_clear (x);
   mpfr_clear (y);
 }
diff --git a/mpfr/tests/tset_str.c b/mpfr/tests/tset_str.c
index cda222d30..5ea2b18d0 100644
--- a/mpfr/tests/tset_str.c
+++ b/mpfr/tests/tset_str.c
@@ -94,6 +94,15 @@ main(int argc, char **argv)
     exit(1);
   }
 
+  mpfr_set_str_raw(x, "+0000");
+  mpfr_set_str_raw(x, "+0000E0");
+  mpfr_set_str_raw(x, "0000E0");
+  if (mpfr_get_d(x) != 0.0) {
+    fprintf(stderr, "Error in mpfr_set_str_raw for s=0.0\n"); 
+    mpfr_clear(x); mpfr_clear(y); 
+    exit(1);
+  }
+
   mpfr_set_str(x, "+243495834958.53452345E1", 10, GMP_RNDN);
   mpfr_set_str(x, "9007199254740993", 10, GMP_RNDN);
   mpfr_set_str(x, "9007199254740992", 10, GMP_RNDU);
diff --git a/mpfr/tests/tsin_cos.c b/mpfr/tests/tsin_cos.c
index 817d10e32..01c7c395f 100644
--- a/mpfr/tests/tsin_cos.c
+++ b/mpfr/tests/tsin_cos.c
@@ -22,8 +22,10 @@ MA 02111-1307, USA. */
 #include "gmp.h"
 #include "mpfr.h"
 
-void large_test (int, int);
-void check53 (double, double, double, mp_rnd_t);
+void large_test _PROTO ((int, int));
+void check53 _PROTO ((double, double, double, mp_rnd_t));
+void check53sin _PROTO ((double, double, mp_rnd_t));
+void check53cos _PROTO ((double, double, mp_rnd_t));
 
 void large_test (int prec, int N)
 {
@@ -70,6 +72,44 @@ void check53 (double x, double sin_x, double cos_x, mp_rnd_t rnd_mode)
   mpfr_clear (c);
 }
 
+void check53sin (double x, double sin_x, mp_rnd_t rnd_mode)
+{
+  mpfr_t xx, s;
+
+  mpfr_init2 (xx, 53);
+  mpfr_init2 (s, 53);
+  mpfr_set_d (xx, x, rnd_mode); /* should be exact */
+  mpfr_sin_cos (s, NULL, xx, rnd_mode);
+  if (mpfr_get_d (s) != sin_x && (!isnan(sin_x) || !isnan(mpfr_get_d(s)))) {
+    fprintf (stderr, "mpfr_sin_cos failed for x=%1.20e, rnd=%s\n", x,
+	     mpfr_print_rnd_mode (rnd_mode));
+    fprintf (stderr, "mpfr_sin_cos gives sin(x)=%1.20e, expected %1.20e\n",
+	     mpfr_get_d (s), sin_x);
+    exit(1);
+  }
+  mpfr_clear (xx);
+  mpfr_clear (s);
+}
+
+void check53cos (double x, double cos_x, mp_rnd_t rnd_mode)
+{
+  mpfr_t xx, c;
+
+  mpfr_init2 (xx, 53);
+  mpfr_init2 (c, 53);
+  mpfr_set_d (xx, x, rnd_mode); /* should be exact */
+  mpfr_sin_cos (NULL, c, xx, rnd_mode);
+  if (mpfr_get_d (c) != cos_x && (!isnan(cos_x) || !isnan(mpfr_get_d(c)))) {
+    fprintf (stderr, "mpfr_sin_cos failed for x=%1.20e, rnd=%s\n", x,
+	     mpfr_print_rnd_mode (rnd_mode));
+    fprintf (stderr, "mpfr_sin_cos gives cos(x)=%1.20e, expected %1.20e\n",
+	     mpfr_get_d (c), cos_x);
+    exit(1);
+  }
+  mpfr_clear (xx);
+  mpfr_clear (c);
+}
+
 /* tsin_cos prec [N] performs N tests with prec bits */
 int main(int argc, char *argv[])
 {
@@ -97,5 +137,10 @@ int main(int argc, char *argv[])
 	   0.537874062022526966409, GMP_RNDZ);
   check53 (1.00591265847407274059,  8.446508805292128885e-1,
 	   0.53531755997839769456,  GMP_RNDN);
+
+  /* check with one argument NULL */
+  check53sin (1.00591265847407274059,  8.446508805292128885e-1, GMP_RNDN);
+  check53cos (1.00591265847407274059, 0.53531755997839769456,  GMP_RNDN);
+
   return 0;
 }
diff --git a/mpfr/trunc.c b/mpfr/trunc.c
index b6f924d35..c5bcf8bdc 100644
--- a/mpfr/trunc.c
+++ b/mpfr/trunc.c
@@ -26,26 +26,27 @@ MA 02111-1307, USA. */
 #include "mpfr.h"
 #include "mpfr-impl.h"
 
-#ifdef OPERATION_floor
+#if MPFR_FLOOR
 #define _MPFR_FLOOR_OR_CEIL
 #define FUNC_NAME mpfr_floor
+#undef MPFR_FLOOR
 #define MPFR_FLOOR 1
 #define MPFR_CEIL 0
 #endif
 
-#ifdef OPERATION_ceil
+#if MPFR_CEIL
 #define _MPFR_FLOOR_OR_CEIL
 #define FUNC_NAME mpfr_ceil
+#undef MPFR_CEIL
 #define MPFR_CEIL 1
 #define MPFR_FLOOR 0
 #endif
 
-#ifdef OPERATION_trunc
+#if MPFR_TRUNC
 #undef FUNC_NAME
 #define FUNC_NAME mpfr_trunc
 #endif
 
-#if 0
 #ifdef _MPFR_FLOOR_OR_CEIL
 static int
 mpn_zero_p (p, n)
@@ -63,7 +64,6 @@ mpn_zero_p (p, n)
   return 1;
 }
 #endif
-#endif
 
 void
 #if __STDC__