* mpfr/*: Update to mpfr-2-0-2-branch 2003-10-18.

11 files changed, 405 insertions, 95 deletions

diff --git a/mpfr/BUGS b/mpfr/BUGS
index 5295f380a..edb69016d 100644
--- a/mpfr/BUGS
+++ b/mpfr/BUGS
@@ -28,11 +28,6 @@ Known bugs:
   partially implemented. For instance, mpfr_pow (z, x, y, rnd) fails for
   very small x and some values of y.
 
-* Elementary functions may return an incorrect ternary value (e.g.
-  in the rounding to nearest mode, when the exact value is close to
-  a representable number). A rewrite of mpfr_can_round (as proposed
-  by Paul Zimmermann) can help to fix this bug.
-
 * The mpfr_set_ld function assumes that the long double type has an
   exponent of at most 15 bits.
 
@@ -73,5 +68,3 @@ Problems due to compiler bugs:
 
   This results in a problem in the mpfr_set_ld function. A workaround is
   to compile set_ld.c with -O0 (no optimization).
-
-  
diff --git a/mpfr/ChangeLog b/mpfr/ChangeLog
index cc55d227e..4e4b56348 100644
--- a/mpfr/ChangeLog
+++ b/mpfr/ChangeLog
@@ -1,5 +1,204 @@
+2003-10-16  Vincent Lefevre  <Vincent.Lefevre@loria.fr>
+
+	* get_d.c, set_d.c: Fixed several bugs.
+
+	* get_d.c: Added XDEBUG support (like in set_d.c).
+
+2003-10-15  Vincent Lefevre  <Vincent.Lefevre@loria.fr>
+
+	* mpfr.texi:
+	Replaced all non-ASCII characters by an equivalent 7-bit sequence
+	to support unpatched texinfo.tex files.
+
+2003-10-15  Paul Zimmermann  <Paul.Zimmermann@loria.fr>
+
+	* tests/reuse.c: removed #define DEBUG (turned on accidentally)
+
+2003-10-15  Vincent Lefevre  <Vincent.Lefevre@loria.fr>
+
+	* BUGS:
+	Bug on the ternary value was fixed by Paul Zimmermann on 2003-10-14.
+
+2003-10-14  Paul Zimmermann  <Paul.Zimmermann@loria.fr>
+
+	* mpfr.texi: added preamble about ternary flag for special functions
+
+	* tests/reuse.c, tests/tacos.c, tests/tagm.c, tests/tasin.c, tests/tatan.c, tests/terf.c, tests/texp.c, tests/tgamma.c, tests/tgeneric.c, tests/thypot.c, tests/tzeta.c:
+	now uses #include "tgeneric.c" everywhere
+	(and modified test_generic to check also the inexact flag)
+
+	* tests/tui_div.c: check_nan() is back
+
+	* tests/tlog.c: removed old unused code
+
+	* tests/thyperbolic.c: fixed tests for x=0
+	removed composition tests (were already in-between #if 0 ... #endif)
+
+	* pow.c, pow_si.c, pow_ui.c, sin.c, sinh.c, sqrt.c, tan.c, tanh.c, ui_pow_ui.c, zeta.c, factorial.c, gamma.c, hypot.c, log.c, log10.c, log1p.c, log2.c, asinh.c, atan.c, atanh.c, const_euler.c, const_log2.c, const_pi.c, cos.c, cosh.c, erf.c, exp2.c, expm1.c, acos.c, acosh.c, agm.c, asin.c:
+	replaced mpfr_can_round (approx, err, rnd1, GMP_RNDN, prec)
+	by       mpfr_can_round (approx, err, rnd1, GMP_RNDZ, prec + 1)
+	which in addition guarantees a correct inexact flag
+
+	* TODO: added new items (version number, rounding modes)
+
+	* README.dev: added comment about --enable-alloca=debug
+
+2003-10-13  Paul Zimmermann  <Paul.Zimmermann@loria.fr>
+
+	* BUGS: added section "Problems due to compiler bugs"
+
+	* exp3.c, exp_2.c:
+	change in can_round calls to get correct inexact flag for rounding to nearest
+
+2003-10-10  Paul Zimmermann  <Paul.Zimmermann@loria.fr>
+
+	* pow.c, tests/tpow.c:
+	fixed bug in mpfr_pow found by Ming J. Tsai (overflow)
+
+2003-10-08  Vincent Lefevre  <Vincent.Lefevre@loria.fr>
+
+	* acos.c, asin.c, atan.c:
+	Removed useless inclusion of standard headers and some #ifdef DEBUG
+	code. Partial reindentation.
+
+	* strncasecmp.c: Added a #include <stddef.h> because size_t is used.
+
+2003-10-07  Paul Zimmermann  <Paul.Zimmermann@loria.fr>
+
+	* tests/tabs.c, tests/tacos.c, tests/tadd.c, tests/tasin.c, tests/tatan.c, tests/tcmp.c, tests/tconst_log2.c, tests/tdiv.c, tests/tfactorial.c, tests/tfrac.c, tests/tgamma.c, tests/tget_str.c, tests/tlog10.c, tests/tlog1p.c, tests/tpow3.c, tests/tset_d.c, tests/tset_si.c, tests/tset_str.c, tests/tset_z.c, tests/tsub.c, tests/tui_pow.c, tests/tzeta.c, tests/tsin.c:
+	reduced test time
+
+2003-10-06  Paul Zimmermann  <Paul.Zimmermann@loria.fr>
+
+	* round_prec.c: fixed comment of mpfr_round_raw_generic
+
+	* add1.c, extract.c, get_si.c, get_ui.c:
+	replaced ABSSIZE by ESIZE (ABSSIZE is the allocated size, and should be
+	used only in functions init, set_prec, round_prec)
+
+	* tests/tadd.c: fixed bug in allocation for in-place operation
+
+	* asin.c, atan.c:
+	inexact flag should now be correct for directed rounding
+	fixed a bug for mpfr_atan(-Inf) [gave +Pi/2 instead of -Pi/2]
+
+	* tests/tatan.c: added test for atan(-Inf)
+
+	* set_ld.c: use macros to avoid possible problem with float input
+
+	* TODO: new proposal for mpfr_can_round
+
+	* mpfr.texi: mpfr_mul_2exp/mpfr_div_2exp are not obsolete
+
+2003-10-05  Vincent Lefevre  <Vincent.Lefevre@loria.fr>
+
+	* mpfr.texi:
+	Clearer mpfr_eq documentation (thanks to Kevin Ryde for the remark).
+
+2003-10-03  Vincent Lefevre  <Vincent.Lefevre@loria.fr>
+
+	* mpfr-impl.h:
+	MPFR_ASSERTN rewritten to avoid "statement with no effect" warnings
+	with gcc when the assertion is always true.
+
+	* add1.c:
+	Optimization (thanks to Patrick Pelissier), as the allocated size
+	may be larger than the size used by the significant bits.
+
+	* mpfr.texi, INSTALL: Updated installation notes.
+
+	* cmp_abs.c, mpfr.texi:
+	Infinities are now accepted in mpfr_cmpabs. Updated its definition
+	in the source (no longer sign(abs(b) - abs(c))).
+
+2003-10-02  Vincent Lefevre  <Vincent.Lefevre@loria.fr>
+
+	* mpfr.texi: Use @url{} for URLs.
+
+	* mpfr.texi: As suggested by PZ:
+	  o documentation of mpfr_fits_* functions moved to the section
+	    on conversions;
+	  o grouped some mpfr_pow special cases;
+	  o explanation concerning the meaning of rnd in mpfr_frac.
+
+	* tests/tset_str.c, tests/tset_z.c, tests/tsinh.c, tests/tsqrt.c, tests/tsub.c, tests/ttanh.c, tests/tui_pow.c, tests/tui_sub.c, tests/tzeta.c, tests/tmul.c, tests/tmul_ui.c, tests/tout_str.c, tests/tpow.c, tests/tset.c, tests/tset_d.c, tests/tlog.c, tests/tlog10.c, tests/tlog1p.c, tests/tlog2.c, tests/texpm1.c, tests/tfma.c, tests/tget_str.c, tests/thypot.c, tests/tadd.c, tests/tadd_ui.c, tests/tagm.c, tests/tasinh.c, tests/tatan.c, tests/tatanh.c, tests/tcan_round.c, tests/tcmp.c, tests/tcmp2.c, tests/tcos.c, tests/tcosh.c, tests/tdiv.c, tests/tdiv_ui.c, tests/terf.c, tests/texceptions.c, tests/texp.c, mpfr.h, mpfr.texi, print_raw.c, set_str_raw.c, tests/tacosh.c, BUGS, mpfr-impl.h:
+	Updated documentation. In particular, mpfr_set_str_raw renamed
+	as mpfr_set_str_binary. This function and mpfr_print_binary are
+	now internal functions. mpfr_print_binary no longer prints the
+	non-significant 0 bits. Updated the source to match the manual.
+	mpfr_print_binary has been completely rewritten (now directly
+	prints to stdout, without using an intermediate string). In
+	mpfr_set_str_binary, replaced atol by strtol + error checking.
+
+2003-10-02  Paul Zimmermann  <Paul.Zimmermann@loria.fr>
+
+	* mpfr.texi: added paragraph on support/grants
+	added help for mpfr_erf
+
+2003-10-02  Vincent Lefevre  <Vincent.Lefevre@loria.fr>
+
+	* mpfr.texi:
+	Make node titles match section titles, and updated menus and xrefs.
+	Added a line break after @samp{uninstall}.
+
+2003-10-01  Vincent Lefevre  <Vincent.Lefevre@loria.fr>
+
+	* TODO, const_log2.c, const_pi.c, mpfr.h, mpfr.texi, zeta.c:
+	Corrections in the MPFR manual (PZ & VL). Functions mpfr_const_pi,
+	mpfr_const_log2 and mpfr_zeta now return a ternary value. Updated
+	TODO file.
+
+2003-09-30  Vincent Lefevre  <Vincent.Lefevre@loria.fr>
+
+	* BUGS: Added a bug related to the ternary value and worst cases.
+
+	* pow_si.c, pow_ui.c: Added exponent range support.
+
+	* mpfr.texi: Corrections up to Section 5.6 (PZ & VL).
+
+	* get_z_exp.c:
+	Re-added the assert on the exponent, but replacing MPFR_EMIN_MIN by
+	MP_EXP_T_MIN (this makes more sense): an assertion failed would mean
+	that the exponent is not representable (an undefined behavior in the
+	ISO C standard). If need be, we could choose to return MP_EXP_T_MIN
+	in such a case, or perhaps MP_EXP_T_MAX to signal an error. The
+	mantissa would still be meaningful.
+
+2003-09-30  Patrick Pelissier  <Patrick.Pelissier@loria.fr>
+
+	* TODO: Update TODO & mpfr_set_prec.
+
+2003-09-29  Vincent Lefevre  <Vincent.Lefevre@loria.fr>
+
+	* mpfr.texi: Corrections up to Section 4.
+
+	* AUTHORS: Added authors Kevin Ryde and Patrick Pelissier.
+
+	* INSTALL: Removed "known problems" that are no longer problems.
+
+2003-09-26  Vincent Lefevre  <Vincent.Lefevre@loria.fr>
+
+	* BUGS: Removed bugs related to the tests and to the exponents
+	as they no longer occur. Updated some potentials bugs.
+
+2003-09-25  Vincent Lefevre  <Vincent.Lefevre@loria.fr>
+
+	* mpfr-impl.h, mpfr.h, mpfr.texi, tests/tinits.c:
+	Removed mpfr_inits, mpfr_inits2, mpfr_clears from the documentation.
+	Moved their prototypes to mpfr-impl.h (internal functions until
+	decided otherwise).
+
+2003-09-25  Patrick Pelissier  <Patrick.Pelissier@loria.fr>
+
+	* const_log2.c, const_pi.c, generic.c, atan.c:
+	Modify 'r' arg of GENERIC from int to long (min 32 bits).
+
 2003-09-25  Vincent Lefevre  <Vincent.Lefevre@loria.fr>
 
+	* mpfr.texi: 8-bit ISO-8859-1 characters + consistent spelling.
+
+	* ChangeLog: Update.
+
 	* tests/tgeneric.c, tests/tget_d.c, tests/tget_str.c, tests/thypot.c, tests/tisnan.c, tests/tlog.c, tests/tlog10.c, tests/tmul.c, tests/tmul_2exp.c, tests/tmul_ui.c, tests/tout_str.c, tests/tpow.c, tests/tpow3.c, tests/trandom.c, tests/tconst_pi.c, tests/tcos.c, tests/tdiv.c, tests/tdiv_ui.c, tests/teq.c, tests/terf.c, tests/tests.c, tests/texceptions.c, tests/texp.c, tests/texp2.c, tests/tfactorial.c, tests/tfma.c, tests/tfrac.c, tests/tgamma.c, tests/mpf_compat.h, tests/reuse.c, tests/tabs.c, tests/tacos.c, tests/tadd.c, tests/tagm.c, tests/tasin.c, tests/tatan.c, tests/tcan_round.c, tests/tcbrt.c, tests/tcmp.c, tests/tcmp2.c, tests/tcmp_d.c, tests/tcmp_ui.c, tests/tconst_log2.c:
 	Changed the remaining stderr to stdout.
 
diff --git a/mpfr/erf.c b/mpfr/erf.c
index 600e0a6d9..8499ab260 100644
--- a/mpfr/erf.c
+++ b/mpfr/erf.c
@@ -29,14 +29,14 @@ MA 02111-1307, USA. */
 
 #define EXP1 2.71828182845904523536 /* exp(1) */
 
-int mpfr_erf_0 _PROTO((mpfr_ptr, mpfr_srcptr, mp_rnd_t)); 
+int mpfr_erf_0 _PROTO((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
 #if 0
-int mpfr_erf_inf _PROTO((mpfr_ptr, mpfr_srcptr, mp_rnd_t)); 
+int mpfr_erf_inf _PROTO((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
 int mpfr_erfc_inf _PROTO((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
 #endif
 
-int 
-mpfr_erf (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) 
+int
+mpfr_erf (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
 {
   double xf;
   int sign_x;
@@ -80,11 +80,15 @@ mpfr_erf (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
 
   if (xf > n * LOG2) /* |erf x| = 1 or 1- */
     {
-      mpfr_set_ui (y, 1, rnd2);
-      if (rnd2 == GMP_RNDD || rnd2 == GMP_RNDZ)
+      if (rnd2 == GMP_RNDN || rnd2 == GMP_RNDU)
         {
-          mpfr_sub_one_ulp (y, GMP_RNDZ); /* 1 - 2^(1-n) */
-          mpfr_add_one_ulp (y, GMP_RNDU);
+          mpfr_set_ui (y, 1, rnd2);
+          inex = 1;
+        }
+      else
+        {
+          mpfr_setmax (y, 0);
+          inex = -1;
         }
     }
   else  /* use Taylor */
@@ -93,9 +97,14 @@ mpfr_erf (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
     }
 
   if (sign_x < 0)
-    MPFR_CHANGE_SIGN(y);
-
-  return inex;
+    {
+      MPFR_CHANGE_SIGN (y);
+      return - inex;
+    }
+  else
+    {
+      return inex;
+    }
 }
 
 /* return x*2^e */
@@ -104,7 +113,7 @@ double mul_2exp (double x, mp_exp_t e)
 {
   if (e > 0)
     {
-      while (e--) 
+      while (e--)
         x *= 2.0;
     }
   else
@@ -123,8 +132,8 @@ double mul_2exp (double x, mp_exp_t e)
    Assumes x is neither NaN nor infinite nor zero.
    Assumes also that e*x^2 <= n (target precision).
  */
-int 
-mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd_mode) 
+int
+mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd_mode)
 {
   mp_prec_t n, m;
   mp_exp_t nuk, sigmak;
@@ -178,7 +187,7 @@ mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd_mode)
             break;
 
           /* tauk <- 1/2 + tauk * 2^sigmak + (1+8k)*2^nuk */
-          tauk = 0.5 + mul_2exp (tauk, sigmak) 
+          tauk = 0.5 + mul_2exp (tauk, sigmak)
             + mul_2exp (1.0 + 8.0 * (double) k, nuk);
         }
 
@@ -224,7 +233,7 @@ mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd_mode)
    f(x) >= f(sqrt(2/e)) ~ 0.7142767512, thus the final partial sum
    should be > 0.5, and MPFR_EXP(s) should always be >= 0.
  */
-int 
+int
 mpfr_erfc_inf (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd)
 {
   mp_prec_t n, m;
@@ -236,7 +245,7 @@ mpfr_erfc_inf (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd)
   double xf = mpfr_get_d1 (x);
 
   n = MPFR_PREC(res); /* target precision */
-  
+
   mpfr_init2 (y, 2);
   mpfr_init2 (s, 2);
   mpfr_init2 (t, 2);
@@ -286,7 +295,7 @@ mpfr_erfc_inf (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd)
             break;
 
           /* tauk <- 1/2 + tauk * 2^sigmak + 2^(2k+2+nuk) */
-          tauk = 0.5 + mul_2exp (tauk, sigmak) 
+          tauk = 0.5 + mul_2exp (tauk, sigmak)
             + mul_2exp (1.0, 2 * k + 2 + nuk);
         }
 
@@ -301,7 +310,7 @@ mpfr_erfc_inf (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd)
       mpfr_sqrt (y, y, GMP_RNDD);
       mpfr_mul (t, t, y, GMP_RNDN);
       mpfr_div (s, s, t, GMP_RNDN);
-      
+
       /* final error bound on s */
       tauk = mul_2exp (3.0, nuk + 5) + 2.0 * tauk + 115.0;
       log2tauk = __gmpfr_ceil_log2 (tauk);
@@ -321,7 +330,7 @@ mpfr_erfc_inf (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd)
    Assumes x is neither NaN nor infinite nor zero.
    Assumes also that e*x^2 > n (target precision).
  */
-int 
+int
 mpfr_erf_inf (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd)
 {
   mp_prec_t n, m;
diff --git a/mpfr/get_d.c b/mpfr/get_d.c
index 42526e46e..85590c0bb 100644
--- a/mpfr/get_d.c
+++ b/mpfr/get_d.c
@@ -28,8 +28,13 @@ MA 02111-1307, USA. */
 #include "mpfr.h"
 #include "mpfr-impl.h"
 
+#ifdef XDEBUG
+#undef _GMP_IEEE_FLOATS
+#endif
 
-static double mpfr_scale2 _PROTO ((double, int));
+#ifndef _GMP_IEEE_FLOATS
+#define _GMP_IEEE_FLOATS 0
+#endif
 
 /* "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.
@@ -40,6 +45,7 @@ static double mpfr_scale2 _PROTO ((double, int));
    compiler+system was seen incorrectly converting from a "float" NaN.  */
 
 #if _GMP_IEEE_FLOATS
+
 /* The "d" field guarantees alignment to a suitable boundary for a double.
    Could use a union instead, if we checked the compiler supports union
    initializers.  */
@@ -67,7 +73,14 @@ 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
+
+#else /* _GMP_IEEE_FLOATS */
+
+#define MPFR_DBL_INFP DBL_POS_INF
+#define MPFR_DBL_INFM DBL_NEG_INF
+#define MPFR_DBL_NAN DBL_NAN
+
+#endif /* _GMP_IEEE_FLOATS */
 
 
 /* multiplies 1/2 <= d <= 1 by 2^exp */
@@ -137,24 +150,12 @@ mpfr_get_d3 (mpfr_srcptr src, mp_exp_t e, mp_rnd_t rnd_mode)
   int negative;
 
   if (MPFR_IS_NAN(src))
-    {
-#ifdef MPFR_DBL_NAN
-      return MPFR_DBL_NAN;
-#else
-      DIVIDE_BY_ZERO;
-#endif
-    }
+    return MPFR_DBL_NAN;
 
   negative = MPFR_SIGN(src) < 0;
 
   if (MPFR_IS_INF(src))
-    {
-#ifdef MPFR_DBL_INFP
-      return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
-#else
-      DIVIDE_BY_ZERO;
-#endif
-    }
+    return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
 
   if (MPFR_IS_ZERO(src))
     return negative ? -0.0 : 0.0;
diff --git a/mpfr/mpfr.texi b/mpfr/mpfr.texi
index 7798ecd90..f7ea1dfa4 100644
--- a/mpfr/mpfr.texi
+++ b/mpfr/mpfr.texi
@@ -181,6 +181,15 @@ log
 @end macro
 @end ifnottex
 
+@c  @pom{} definition
+@tex
+\gdef\pom{\ifmmode\pm\else$\pm$\fi}
+@end tex
+@ifnottex
+@macro pom
+±
+@end macro
+@end ifnottex
 
 @node Copying, Introduction to MPFR, Top, Top
 @comment  node-name, next, previous,  up
@@ -1026,11 +1035,11 @@ rounded in the direction @var{rnd}.
 Special values are currently handled as described in the ISO C99 standard
 for the @code{pow} function:
 @itemize @bullet
-@item @code{pow(±0, @var{y})} returns plus or minus infinity for @var{y} a negative odd integer.
-@item @code{pow(±0, @var{y})} returns plus infinity for @var{y} negative and not an odd integer.
-@item @code{pow(±0, @var{y})} returns plus or minus zero for @var{y} a positive odd integer.
-@item @code{pow(±0, @var{y})} returns plus zero for @var{y} positive and not an odd integer.
-@item @code{pow(-1, ±inf)} returns 1.
+@item @code{pow(@pom{}0, @var{y})} returns plus or minus infinity for @var{y} a negative odd integer.
+@item @code{pow(@pom{}0, @var{y})} returns plus infinity for @var{y} negative and not an odd integer.
+@item @code{pow(@pom{}0, @var{y})} returns plus or minus zero for @var{y} a positive odd integer.
+@item @code{pow(@pom{}0, @var{y})} returns plus zero for @var{y} positive and not an odd integer.
+@item @code{pow(-1, @pom{}inf)} returns 1.
 @item @code{pow(+1, @var{y})} returns 1 for any @var{x}, even a NaN.
 @item @code{pow(@var{x}, @var{y})} returns NaN for finite negative @var{x} and finite non-integer @var{y}.
 @item @code{pow(@var{x}, -inf)} returns plus infinity for @math{0 < @GMPabs{x} < 1}, and plus zero for @math{@GMPabs{x} > 1}.
@@ -1537,11 +1546,11 @@ in raw binary format (the exponent is written in decimal, yet).
 @comment  node-name,  next,  previous,  up
 @unnumbered Contributors
 
-The main developers consist of Guillaume Hanrot, Vincent Lefèvre,
+The main developers consist of Guillaume Hanrot, Vincent Lef@`evre,
 Kevin Ryde and Paul Zimmermann.
 
 We would like to thank Jean-Michel Muller and Joris van der Hoeven for very
-fruitful discussions at the beginning of that project, Torbjörn Granlund
+fruitful discussions at the beginning of that project, Torbj@"orn Granlund
 and Kevin Ryde
 for their help about design issues
 and their suggestions for an easy integration into GNU MP,
@@ -1561,13 +1570,13 @@ David Daney contributed the hyperbolic and inverse hyperbolic functions,
 the base-2 exponential, and the factorial function. Fabrice Rouillier
 contributed the original version of @file{mul_ui.c}, the @file{gmp_op.c}
 file, and helped to the Windows porting.
-Jean-Luc Rémy contributed the @code{mpfr_zeta} code.
+Jean-Luc R@'emy contributed the @code{mpfr_zeta} code.
 Ludovic Meunier helped in the design of the @code{mpfr_erf} code.
 
 The development of the MPFR library would not have been possible without the
 continuous support of LORIA, INRIA and INRIA Lorraine.
 The development of MPFR was also supported by a grant 
-(202F0659 00 MPN 121) from the Conseil Régional de Lorraine in 2002.
+(202F0659 00 MPN 121) from the Conseil R@'egional de Lorraine in 2002.
 
 @node References, GNU Free Documentation License, Contributors, Top
 @comment  node-name,  next,  previous,  up
@@ -1576,7 +1585,7 @@ The development of MPFR was also supported by a grant
 @itemize @bullet
 
 @item
-Torbjörn Granlund, "GNU MP: The GNU Multiple Precision Arithmetic Library",
+Torbj@"orn Granlund, "GNU MP: The GNU Multiple Precision Arithmetic Library",
   version 4.1.2, 2002.
 
 @item
diff --git a/mpfr/set_d.c b/mpfr/set_d.c
index a98256d3c..a6fa6bf42 100644
--- a/mpfr/set_d.c
+++ b/mpfr/set_d.c
@@ -30,12 +30,10 @@ MA 02111-1307, USA. */
 #define MPFR_LIMBS_PER_DOUBLE 2
 #elif (BITS_PER_MP_LIMB >= 64)
 #define MPFR_LIMBS_PER_DOUBLE 1
-#elif (BITS_PER_MP_LIMB == 16)
-#define MPFR_LIMBS_PER_DOUBLE 4
+#else
+#error "Unsupported value of BITS_PER_MP_LIMB"
 #endif
 
-static int __mpfr_extract_double _PROTO ((mp_ptr, double));
-
 /* Included from gmp-2.0.2, patched to support denorms */
 
 #ifdef XDEBUG
@@ -57,10 +55,9 @@ __mpfr_extract_double (mp_ptr rp, double d)
 #endif
 
   /* BUGS
-
      1. Should handle Inf and NaN in IEEE specific code.
      2. Handle Inf and NaN also in default code, to avoid hangs.
-     3. Generalize to handle all BITS_PER_MP_LIMB >= 32.
+     3. Generalize to handle all BITS_PER_MP_LIMB.
      4. This lits is incomplete and misspelled.
    */
 
@@ -71,6 +68,7 @@ __mpfr_extract_double (mp_ptr rp, double d)
     }
 
 #if _GMP_IEEE_FLOATS
+
   {
     union ieee_double_extract x;
     x.d = d;
@@ -96,16 +94,21 @@ __mpfr_extract_double (mp_ptr rp, double d)
 	manl = x.s.manl << 11; /* low 21 bits */
 #endif
       }
+
+    if (exp)
+      exp -= 1022;
+    else
+      exp = -1021;
   }
-#else
+
+#else /* _GMP_IEEE_FLOATS */
+
   {
     /* Unknown (or known to be non-IEEE) double format.  */
     exp = 0;
     if (d >= 1.0)
       {
-        if (d * 0.5 == d)
-          abort ();
-
+        MPFR_ASSERTN (d * 0.5 != d);
         while (d >= 32768.0)
           {
             d *= (1.0 / 65536.0);
@@ -138,12 +141,9 @@ __mpfr_extract_double (mp_ptr rp, double d)
     manh = d;
     manl = (d - manh) * MP_BASE_AS_DOUBLE;
 #endif
-
-    exp += 1022;
   }
-#endif
 
-  if (exp) exp = (unsigned) exp - 1022; else exp = -1021;
+#endif /* _GMP_IEEE_FLOATS */
 
 #if BITS_PER_MP_LIMB == 64
   rp[0] = manl;
diff --git a/mpfr/tests/reuse.c b/mpfr/tests/reuse.c
index b410b94a3..e11fe3d40 100644
--- a/mpfr/tests/reuse.c
+++ b/mpfr/tests/reuse.c
@@ -46,8 +46,6 @@ mpfr_compare (mpfr_t a, mpfr_t b)
     (MPFR_IS_NAN(b) || mpfr_cmp(a, b));
 }
 
-#define DEBUG
-
 void
 test3 (char *foo, mp_prec_t prec, mp_rnd_t rnd)
 {
diff --git a/mpfr/tests/tasin.c b/mpfr/tests/tasin.c
index 2aef35121..bba0c5221 100644
--- a/mpfr/tests/tasin.c
+++ b/mpfr/tests/tasin.c
@@ -23,7 +23,9 @@ MA 02111-1307, USA. */
 #include <stdio.h>
 #include <stdlib.h>
 #include "gmp.h"
+#include "gmp-impl.h"
 #include "mpfr.h"
+#include "mpfr-impl.h"
 #include "mpfr-test.h"
 
 #define TEST_FUNCTION mpfr_asin
diff --git a/mpfr/tests/terf.c b/mpfr/tests/terf.c
index a159864ec..64c1cca72 100644
--- a/mpfr/tests/terf.c
+++ b/mpfr/tests/terf.c
@@ -36,6 +36,7 @@ int
 main (int argc, char *argv[])
 {
   mpfr_t x, y;
+  int inex;
 
   tests_start_mpfr ();
 
@@ -159,6 +160,44 @@ main (int argc, char *argv[])
       exit (1);
     }
 
+  mpfr_set_prec (x, 8);
+  mpfr_set_prec (y, 8);
+  mpfr_set_ui (x, 50, GMP_RNDN);
+  inex = mpfr_erf (y, x, GMP_RNDN);
+  if (mpfr_cmp_ui (y, 1))
+    {
+      printf ("mpfr_erf failed for x=50, rnd=GMP_RNDN\n");
+      printf ("expected 1, got ");
+      mpfr_out_str (stdout, 2, 0, y, GMP_RNDN);
+      printf ("\n");
+      exit (1);
+    }
+  if (inex <= 0)
+    {
+      printf ("mpfr_erf failed for x=50, rnd=GMP_RNDN: wrong ternary value\n"
+              "expected positive, got %d\n", inex);
+      exit (1);
+    }
+  inex = mpfr_erf (x, x, GMP_RNDZ);
+  mpfr_nextbelow (y);
+  if (mpfr_cmp (x, y))
+    {
+      printf ("mpfr_erf failed for x=50, rnd=GMP_RNDZ\n");
+      printf ("expected ");
+      mpfr_out_str (stdout, 2, 0, y, GMP_RNDN);
+      printf ("\n");
+      printf ("got      ");
+      mpfr_out_str (stdout, 2, 0, x, GMP_RNDN);
+      printf ("\n");
+      exit (1);
+    }
+  if (inex >= 0)
+    {
+      printf ("mpfr_erf failed for x=50, rnd=GMP_RNDN: wrong ternary value\n"
+              "expected negative, got %d\n", inex);
+      exit (1);
+    }
+
   mpfr_clear (x);
   mpfr_clear (y);
 
diff --git a/mpfr/tests/tgamma.c b/mpfr/tests/tgamma.c
index c3179e32b..b19621583 100644
--- a/mpfr/tests/tgamma.c
+++ b/mpfr/tests/tgamma.c
@@ -22,7 +22,9 @@ MA 02111-1307, USA. */
 #include <stdio.h>
 #include <stdlib.h>
 #include "gmp.h"
+#include "gmp-impl.h"
 #include "mpfr.h"
+#include "mpfr-impl.h"
 #include "mpfr-test.h"
 
 int mpfr_gamma (mpfr_ptr, mpfr_srcptr, mp_rnd_t);
diff --git a/mpfr/tests/trint.c b/mpfr/tests/trint.c
index e8a47cefa..7c4d44d3f 100644
--- a/mpfr/tests/trint.c
+++ b/mpfr/tests/trint.c
@@ -1,4 +1,4 @@
-/* Test file for mpfr_trunc, mpfr_floor, mpfr_ceil, mpfr_round.
+/* Test file for mpfr_rint, mpfr_trunc, mpfr_floor, mpfr_ceil, mpfr_round.
 
 Copyright 2002, 2003 Free Software Foundation.
 
@@ -25,6 +25,56 @@ MA 02111-1307, USA. */
 #include "mpfr.h"
 #include "mpfr-test.h"
 
+#if __STDC_VERSION__ >= 199901L
+
+static void
+test_fct (double (*f)(double), int (*g)(), char *s, mp_rnd_t r)
+{
+  double d, y;
+  mpfr_t dd, yy;
+
+  mpfr_init2 (dd, 53);
+  mpfr_init2 (yy, 53);
+  for (d = -5.0; d <= 5.0; d += 0.25)
+    {
+      mpfr_set_d (dd, d, r);
+      y = (*f)(d);
+      if (g == &mpfr_rint)
+        mpfr_rint (yy, dd, r);
+      else
+        (*g)(yy, dd);
+      mpfr_set_d (dd, y, r);
+      if (mpfr_cmp (yy, dd))
+        {
+          printf ("test_against_libc: incorrect result for %s, rnd = %s,"
+                  " d = %g\ngot ", s, mpfr_print_rnd_mode (r), d);
+          mpfr_out_str (stdout, 10, 0, yy, GMP_RNDN);
+          printf (" instead of %g\n", y);
+          exit (1);
+        }
+    }
+  mpfr_clear (dd);
+  mpfr_clear (yy);
+}
+
+#define TEST_FCT(F) test_fct (&F, &mpfr_##F, #F, r)
+
+static void
+test_against_libc (void)
+{
+  int r = 0;
+
+  TEST_FCT (round);
+  TEST_FCT (trunc);
+  TEST_FCT (floor);
+  TEST_FCT (ceil);
+  for (r = 0; r < 4; r++)
+    if (mpfr_set_machine_rnd_mode (r) == 0)
+      TEST_FCT (rint);
+}
+
+#endif
+
 int
 main (void)
 {
@@ -57,33 +107,37 @@ main (void)
                   (unsigned int) s);
           exit (1);
         }
-      for (p=2; p<100; p++)
+      for (p = 2; p < 100; p++)
         {
+          int trint;
           mpfr_set_prec (y, p);
-          for (r=0; r<4; r++)
-            {
-              if (r == GMP_RNDN)
-                inexact = mpfr_round (y, x);
-              else if (r == GMP_RNDZ)
-                inexact = mpfr_trunc (y, x);
-              else if (r == GMP_RNDU)
-                inexact = mpfr_ceil (y, x);
-              else /* r = GMP_RNDD */
-                inexact = mpfr_floor (y, x);
-              if (mpfr_sub (t, y, x, GMP_RNDN))
-                {
-                  printf ("Error: subtraction should be exact\n");
-                  exit (1);
-                }
-              sign_t = mpfr_cmp_ui (t, 0);
-              if (((inexact == 0) && (sign_t != 0)) ||
-                  ((inexact < 0) && (sign_t >= 0)) ||
-                  ((inexact > 0) && (sign_t <= 0)))
-                {
-                  printf ("Wrong inexact flag\n");
-                  exit (1);
-                }
-            }
+          for (r = 0; r < 4; r++)
+            for (trint = 0; trint < 2; trint++)
+              {
+                if (trint)
+                  inexact = mpfr_rint (y, x, r);
+                else if (r == GMP_RNDN)
+                  inexact = mpfr_round (y, x);
+                else if (r == GMP_RNDZ)
+                  inexact = mpfr_trunc (y, x);
+                else if (r == GMP_RNDU)
+                  inexact = mpfr_ceil (y, x);
+                else /* r = GMP_RNDD */
+                  inexact = mpfr_floor (y, x);
+                if (mpfr_sub (t, y, x, GMP_RNDN))
+                  {
+                    printf ("Error: subtraction should be exact\n");
+                    exit (1);
+                  }
+                sign_t = mpfr_cmp_ui (t, 0);
+                if (((inexact == 0) && (sign_t != 0)) ||
+                    ((inexact < 0) && (sign_t >= 0)) ||
+                    ((inexact > 0) && (sign_t <= 0)))
+                  {
+                    printf ("Wrong inexact flag\n");
+                    exit (1);
+                  }
+              }
         }
     }
   mpfr_clear (x);
@@ -91,6 +145,10 @@ main (void)
   mpz_clear (z);
   mpfr_clear (t);
 
+#if __STDC_VERSION__ >= 199901L
+  test_against_libc ();
+#endif
+
   tests_end_mpfr ();
   return 0;
 }