added mpfr_remquo and mpfr_remainder

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@4420 280ebfd0-de03-0410-8827-d642c229c3f4

7 files changed, 438 insertions, 15 deletions

diff --git a/Makefile.am b/Makefile.am
index 13ab1a4ba..f5cf9cd2c 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -7,7 +7,7 @@ include_HEADERS = mpfr.h mpf2mpfr.h
 
 lib_LTLIBRARIES = libmpfr.la
 
-libmpfr_la_SOURCES = mpfr.h mpf2mpfr.h mpfr-gmp.h mpfr-impl.h mpfr-longlong.h mpfr-test.h exceptions.c extract.c uceil_exp2.c uceil_log2.c ufloor_log2.c add.c add1.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_si.c cmp_ui.c comparisons.c div_2exp.c div_2si.c div_2ui.c div.c div_ui.c dump.c eq.c exp10.c exp2.c exp3.c exp.c frac.c get_d.c get_exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c const_log2.c log.c mul_2exp.c mul_2si.c mul_2ui.c mul.c mul_ui.c neg.c next.c out_str.c const_pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c round_prec.c set.c setmax.c setmin.c set_d.c set_dfl_prec.c set_exp.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c set_str.c set_str_raw.c set_ui.c set_z.c sqrt.c sqrt_ui.c sub.c sub1.c sub_ui.c rint.c ui_div.c ui_sub.c urandomb.c get_z_exp.c swap.c  factorial.c cosh.c sinh.c tanh.c acosh.c asinh.c atanh.c atan.c cmp2.c exp_2.c asin.c const_euler.c cos.c sin.c tan.c fma.c hypot.c log1p.c expm1.c log2.c log10.c ui_pow.c ui_pow_ui.c minmax.c dim.c copysign.c gmp_op.c init2.c acos.c sin_cos.c set_nan.c set_inf.c powerof2.c gamma.c set_ld.c get_ld.c cbrt.c volatile.c fits_s.h fits_sshort.c fits_sint.c fits_slong.c fits_u.h fits_ushort.c fits_uint.c fits_ulong.c fits_uintmax.c fits_intmax.c get_si.c get_ui.c zeta.c cmp_d.c erf.c inits.c inits2.c clears.c sgn.c check.c sub1sp.c version.c mpn_exp.c mpfr-gmp.c mp_clz_tab.c sum.c add1sp.c free_cache.c si_op.c cmp_ld.c set_ui_2exp.c set_si_2exp.c set_uj.c set_sj.c get_sj.c get_uj.c get_z.c iszero.c cache.c sqr.c int_ceil_log2.c isqrt.c strtofr.c pow_z.c logging.c mulders.c get_f.c round_p.c erfc.c atan2.c subnormal.c const_catalan.c root.c gen_inverse.h sec.c csc.c cot.c eint.c sech.c csch.c coth.c round_near_x.c constant.c abort_prec_max.c stack_interface.c lngamma.c zeta_ui.c set_d64.c get_d64.c jn.c yn.c
+libmpfr_la_SOURCES = mpfr.h mpf2mpfr.h mpfr-gmp.h mpfr-impl.h mpfr-longlong.h mpfr-test.h exceptions.c extract.c uceil_exp2.c uceil_log2.c ufloor_log2.c add.c add1.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_si.c cmp_ui.c comparisons.c div_2exp.c div_2si.c div_2ui.c div.c div_ui.c dump.c eq.c exp10.c exp2.c exp3.c exp.c frac.c get_d.c get_exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c const_log2.c log.c mul_2exp.c mul_2si.c mul_2ui.c mul.c mul_ui.c neg.c next.c out_str.c const_pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c round_prec.c set.c setmax.c setmin.c set_d.c set_dfl_prec.c set_exp.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c set_str.c set_str_raw.c set_ui.c set_z.c sqrt.c sqrt_ui.c sub.c sub1.c sub_ui.c rint.c ui_div.c ui_sub.c urandomb.c get_z_exp.c swap.c  factorial.c cosh.c sinh.c tanh.c acosh.c asinh.c atanh.c atan.c cmp2.c exp_2.c asin.c const_euler.c cos.c sin.c tan.c fma.c hypot.c log1p.c expm1.c log2.c log10.c ui_pow.c ui_pow_ui.c minmax.c dim.c copysign.c gmp_op.c init2.c acos.c sin_cos.c set_nan.c set_inf.c powerof2.c gamma.c set_ld.c get_ld.c cbrt.c volatile.c fits_s.h fits_sshort.c fits_sint.c fits_slong.c fits_u.h fits_ushort.c fits_uint.c fits_ulong.c fits_uintmax.c fits_intmax.c get_si.c get_ui.c zeta.c cmp_d.c erf.c inits.c inits2.c clears.c sgn.c check.c sub1sp.c version.c mpn_exp.c mpfr-gmp.c mp_clz_tab.c sum.c add1sp.c free_cache.c si_op.c cmp_ld.c set_ui_2exp.c set_si_2exp.c set_uj.c set_sj.c get_sj.c get_uj.c get_z.c iszero.c cache.c sqr.c int_ceil_log2.c isqrt.c strtofr.c pow_z.c logging.c mulders.c get_f.c round_p.c erfc.c atan2.c subnormal.c const_catalan.c root.c gen_inverse.h sec.c csc.c cot.c eint.c sech.c csch.c coth.c round_near_x.c constant.c abort_prec_max.c stack_interface.c lngamma.c zeta_ui.c set_d64.c get_d64.c jn.c yn.c remquo.c
 
 libmpfr_la_LIBADD = @LIBOBJS@
 
diff --git a/TODO b/TODO
index 6e7d34fb0..9c8c83296 100644
--- a/TODO
+++ b/TODO
@@ -110,19 +110,6 @@ New functions to implement:
      mpfr_clear (q);
      return inexact;
   }
-- mpfr_remainder (mpfr_t, mpfr_srcptr, mpfr_srcptr, mp_rnd_t)
-  as in ISO C99: "centered" remainder.
-- mpfr_remquo (mpfr_t, mpfr_t, mpfr_srcptr, mpfr_srcptr, mp_rnd_t)
-  as in ISO C99: remainder and quotient [suggested by Kaveh Ghazi, 18 Feb 2007]
-  Question: compute the full quotient (mpfr or mpz?) or only its low
-  significant bits as in C99 (see rationale below)?
-  Vincent suggests to store the low bits of the quotient in an mpfr_t,
-  whose precision would be chosen by the user.
-    "New feature for C99. The remquo functions are intended for
-    implementing argument reductions which can exploit a few low-order
-    bits of the quotient. Note that x may be so large in magnitude
-    relative to y that an exact representation of the quotient is not
-    practical."
 - mpfr_fms for a-b*c
   [suggested by Tomas Zahradnicky <tomas@24uSoftware.com>, 29 Nov 2003]
 - 1/sqrt(x) [Regis Dupont <dupont@lix.polytechnique.fr>, 15 Sep 2004]
diff --git a/mpfr.h b/mpfr.h
index cb607a3c9..05838204a 100644
--- a/mpfr.h
+++ b/mpfr.h
@@ -446,6 +446,10 @@ __MPFR_DECLSPEC int mpfr_rint_ceil _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
 __MPFR_DECLSPEC int mpfr_rint_floor _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
                                                   mpfr_rnd_t));
 __MPFR_DECLSPEC int mpfr_frac _MPFR_PROTO ((mpfr_ptr,mpfr_srcptr,mpfr_rnd_t));
+__MPFR_DECLSPEC int mpfr_remquo _MPFR_PROTO ((int*, mpfr_ptr, mpfr_srcptr,
+                                              mpfr_srcptr, mp_rnd_t));
+__MPFR_DECLSPEC int mpfr_remainder _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
+                                                 mpfr_srcptr, mp_rnd_t));
 
 __MPFR_DECLSPEC int mpfr_fits_ulong_p _MPFR_PROTO((mpfr_srcptr, mpfr_rnd_t));
 __MPFR_DECLSPEC int mpfr_fits_slong_p _MPFR_PROTO((mpfr_srcptr, mpfr_rnd_t));
diff --git a/mpfr.texi b/mpfr.texi
index c9dd9a583..9f548d85d 100644
--- a/mpfr.texi
+++ b/mpfr.texi
@@ -1687,6 +1687,19 @@ Set @var{rop} to the fractional part of @var{op}, having the same sign as
 the fractional part is generated).
 @end deftypefun
 
+@deftypefun int mpfr_remainder (mpfr_t @var{r}, mpfr_t @var{x}, mpfr_t @var{y}, mp_rnd_t @var{rnd})
+@deftypefunx int mpfr_remquo (int* @var{q}, mpfr_t @var{r}, mpfr_t @var{x}, mpfr_t @var{y}, mp_rnd_t @var{rnd})
+Set @var{r} to the remainder of the division of @var{x} by @var{y}, with 
+quotient rounded to the nearest integer (ties rounded to even), and
+@var{r} rounded according to the direction @var{rnd}.
+The return value is the inexact flag corresponding to @var{r}.
+Additionally, @code{mpfr_remquo} stores
+the 3 low significant bits from the quotient in @var{q[0]}.
+Note that @var{x} may be so large in magnitude relative to @var{y} that an
+exact representation of the quotient is not practical.
+These functions are useful for additive argument reduction.
+@end deftypefun
+
 @deftypefun int mpfr_integer_p (mpfr_t @var{op})
 Return non-zero iff @var{op} is an integer.
 @end deftypefun
diff --git a/remquo.c b/remquo.c
new file mode 100644
index 000000000..3405cc9ae
--- /dev/null
+++ b/remquo.c
@@ -0,0 +1,208 @@
+/* mpfr_remquo and mpfr_remainder -- argument reduction functions
+
+Copyright 2007 Free Software Foundation, Inc.
+Contributed by the Arenaire and Cacao projects, INRIA.
+
+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., 51 Franklin St, Fifth Floor, Boston,
+MA 02110-1301, USA. */
+
+#include "mpfr-impl.h"
+
+/*
+  Let q = x/y rounded to the nearest integer (to the nearest even number
+  in case x/y = n + 1/2 with n integer).
+  Put x - q*y in rem, rounded according to rnd.
+  The value stored in *quo has the sign of q, and agrees with q with (at least)
+  the 3 low order bits. In other words, *quo = q (mod 8) and *quo q >= 0.
+  If rem is zero, then it has the sign of x.
+  The returned 'int' is the inexact flag giving the place of rem wrt x - q*y.
+
+  If x or y is NaN: *quo is undefined, rem is NaN.
+  If x is Inf, whatever y: *quo is undefined, rem is NaN.
+  If y is Inf, x not NaN nor Inf: *quo is 0, rem is x.
+  If y is 0, whatever x: *quo is undefined, rem is NaN.
+  If x is 0, whatever y (not NaN nor 0): *quo is 0, rem is x.
+
+  Otherwise if x and y are neither NaN, Inf nor 0, q is always defined,
+  thus *quo is.
+  Since |x - q*y| <= y/2, no overflow is possible.
+  Only an underflow is possible when y is very small.
+ */
+
+/* compares the (absolute values of the) significands of x and y,
+   returns a positive value if m(x) > m(y),
+           zero if m(x) = m(y),
+           a negative value of m(x) < m(y),
+   where 1/2 <= m(x), m(y) < 1.
+   Assumes x and y are neither NaN, nor Inf, nor zero. */
+static int
+mpfr_cmp_significand (mpfr_srcptr x, mpfr_srcptr y)
+{
+  mp_prec_t xn = MPFR_LIMB_SIZE (x), yn = MPFR_LIMB_SIZE (y);
+  mp_ptr    xp = MPFR_MANT (x),      yp = MPFR_MANT (y);
+  int i;
+
+  if (xn <= yn)
+    {
+      i = mpn_cmp (xp, yp + (yn - xn), xn);
+      if (i != 0)
+        return i;
+      /* if i = 0, either m(x) = m(y) if {yp, yn - xn} is zero,
+         or m(x) < m(y) */
+      for (yn = yn - xn; yn > 0 && yp[yn - 1] == 0; yn --);
+      return (yn == 0) ? 0 : -1;
+    }
+  else /* xn > yn */
+    {
+      i = mpn_cmp (xp + (xn - yn), yp, yn);
+      if (i != 0)
+        return i;
+      /* i = 0: need to check the low xn-yn limbs from x */
+      for (xn = xn - yn; xn > 0 && xp[xn - 1] == 0; xn--);
+      /* if xn > 0, one of the low limbs of x is non-zero */
+      return (xn == 0) ? 0 : 1;
+    }
+}
+
+#define WANTED_BITS 3
+
+#if (WANTED_BITS < 3)
+#error "WANTED_BITS must be at least 3 to conform to C99"
+#endif
+
+#if (WANTED_BITS > BITS_PER_MP_LIMB)
+#error "WANTED_BITS must be less or equal to BITS_PER_MP_LIMB"
+#endif
+
+/* assuming q is an integer, with in addition EXP(q) <= PREC(q),
+   returns |q| mod 2^WANTED_BITS */
+static int
+get_low_bits (mpfr_srcptr q)
+{
+  mp_ptr qp = MPFR_MANT(q);
+  mp_size_t qn = MPFR_LIMB_SIZE(q);
+  mp_size_t w = qn * BITS_PER_MP_LIMB - MPFR_EXP(q);
+  mp_limb_t res;
+
+  /* weight of bit 0 of q is -w, with -w <= 0 since EXP(q) <= PREC(q),
+     thus bit of weight 0 is w. */
+  /* normally this loop should not be used, since in normal cases we have
+     ulp(q)=1, i.e., EXP(q) = PREC(q), thus w is exactly the number of
+     unused bits in the low significant limb of q. */
+  while (w >= BITS_PER_MP_LIMB)
+    {
+      qp ++;
+      qn --;
+      w -= BITS_PER_MP_LIMB;
+    }
+  if ((w + WANTED_BITS <= BITS_PER_MP_LIMB) || (qn <= 1))
+    /* all wanted bits in same limb */
+    res = qp[0] >> w;
+  else /* wanted bits are split in two consecutive limbs: BITS_PER_MP_LIMB-w
+          in qp[0], and WANTED_BITS-(BITS_PER_MP_LIMB-w) in qp[1] */
+    res = (qp[0] >> w) | (qp[1] << (BITS_PER_MP_LIMB - w));
+  return res & ((MPFR_LIMB_ONE << WANTED_BITS) - MPFR_LIMB_ONE);
+}
+
+int
+mpfr_remquo (int *quo, mpfr_ptr rem,
+             mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd)
+{
+  mpfr_t q, qy;
+  int inex, compare;
+  mp_exp_t ex, ey;
+  mp_prec_t prec_q;
+
+  MPFR_LOG_FUNC (("x[%#R]=%R y[%#R]=%R rnd=%d", x, x, y, y, r),
+                 ("quo=%d rem[%#R]=%R", *quo, rem, rem));
+
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x) || MPFR_IS_SINGULAR (y)))
+    {
+      if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y) || MPFR_IS_INF (x)
+          || MPFR_IS_ZERO (y))
+        {
+          MPFR_SET_NAN (rem);
+          MPFR_RET_NAN;
+        }
+      else /* either y is Inf and x is 0 or non-special,
+              or x is 0 and y is non-special,
+              in both cases the quotient is zero. */
+        {
+          *quo = 0;
+          return mpfr_set (rem, x, rnd);
+        }
+    }
+
+  /* now neither x nor y is NaN, Inf or zero */
+  
+  /* first deal with the easy case where x is already reduced mod y,
+     i.e., |x| <= |y|/2 */
+  ex = MPFR_EXP (x);
+  ey = MPFR_EXP (y);
+  compare = mpfr_cmp_significand (x, y); /* compare > 0 if m(x) > m(y) */
+  if ((ex + 1 < ey) || ((ex + 1 == ey) && compare <= 0))
+    {
+      *quo = 0;
+      return mpfr_set (rem, x, rnd);
+    }
+
+  /* Now we are sure that the (true) quotient q is > 1/2 in absolute value.
+     If ex = EXP(x) and ey = EXP(y), we have |x| < 2^ex and |y| >= 2^(ey-1),
+     thus |q| < 2^(ex-ey+1).
+     We also know that ex + 1 >= ey. */
+
+  prec_q = ex - ey + (compare >= 0);
+  mpfr_init2 (q, (prec_q < MPFR_PREC_MIN) ? MPFR_PREC_MIN : prec_q);
+  inex = mpfr_div (q, x, y, GMP_RNDN);
+  if (prec_q < MPFR_PREC_MIN)
+    {
+      int inex2;
+      /* we might have a double-rounding problem in case q = n + 1/2 here,
+         which can be obtained either from x/y = n + 1/2 + eps, or
+         x/y = n + 1/2 - eps */
+      inex2 = mpfr_rint (q, q, GMP_RNDN);
+      if (inex2 == MPFR_EVEN_INEX && inex > 0)
+        /* mpfr_div rounded n + 1/2 - eps to n + 1/2, and mpfr_rint
+           rounded n + 1/2 to n + 1 */
+        mpfr_sub_ui (q, q, 1, GMP_RNDN);
+      else if (inex2 == -MPFR_EVEN_INEX && inex < 0)
+        mpfr_add_ui (q, q, 1, GMP_RNDN);
+    }
+
+  /* set the low bits of the quotient */
+  *quo = get_low_bits (q);
+  /* quotient should have the sign of x/y */
+  if (MPFR_SIGN(x) != MPFR_SIGN(y))
+    *quo = -*quo;
+
+  /* since we have no fused-multiply-and-subtract yet, we compute
+     x + (-q)*y with an FMA */
+  mpfr_neg (q, q, GMP_RNDN); /* exact */
+  inex = mpfr_fma (rem, q, y, x, rnd);
+
+  mpfr_clear (q);
+
+  return inex;
+}
+
+int
+mpfr_remainder (mpfr_ptr rem, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd)
+{
+  int quo[1];
+
+  return mpfr_remquo (quo, rem, x, y, rnd);
+}
diff --git a/tests/Makefile.am b/tests/Makefile.am
index dee71e676..6549026d4 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -1,6 +1,6 @@
 AUTOMAKE_OPTIONS = 1.6 gnu $(top_builddir)/ansi2knr
 
-check_PROGRAMS = tversion tinternals tinits tsgn tcheck tisnan texceptions tset_exp tset tabs tset_d tset_f tset_q tset_si tset_str tset_z tset_ld tset_sj tswap tcopysign tcmp tcmp2 tcmpabs tcmp_d tcmp_ld tcomparisons teq tadd tsub tmul tdiv tsub1sp tadd1sp tadd_ui tsub_ui tcmp_ui tdiv_ui tmul_ui tsqrt_ui tui_div tui_sub tgmpop tsi_op tmul_2exp tfma tsum tdim tminmax tnext tfits tget_d tget_d_2exp tget_z tget_str tget_sj tout_str tinp_str toutimpl tcan_round tround_prec tsqrt tconst_log2 tconst_pi tconst_euler trandom ttrunc trint tfrac texp texp2 texpm1 tlog tlog2 tlog10 tlog1p tpow tui_pow tpow3 tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic tasin tacos tcos tatan tsin ttan tsin_cos tagm thypot tfactorial tgamma terf tcbrt tzeta mpf_compat mpfr_compat reuse tsqr tstrtofr tpow_z tget_f tconst_catalan troot tsec tcsc tcot teint tcoth tcsch tsech tstckintc tsubnormal tlngamma tlgamma tzeta_ui tget_ld_2exp tget_set_d64 tj0 tj1 tjn ty0 ty1 tyn
+check_PROGRAMS = tversion tinternals tinits tsgn tcheck tisnan texceptions tset_exp tset tabs tset_d tset_f tset_q tset_si tset_str tset_z tset_ld tset_sj tswap tcopysign tcmp tcmp2 tcmpabs tcmp_d tcmp_ld tcomparisons teq tadd tsub tmul tdiv tsub1sp tadd1sp tadd_ui tsub_ui tcmp_ui tdiv_ui tmul_ui tsqrt_ui tui_div tui_sub tgmpop tsi_op tmul_2exp tfma tsum tdim tminmax tnext tfits tget_d tget_d_2exp tget_z tget_str tget_sj tout_str tinp_str toutimpl tcan_round tround_prec tsqrt tconst_log2 tconst_pi tconst_euler trandom ttrunc trint tfrac texp texp2 texpm1 tlog tlog2 tlog10 tlog1p tpow tui_pow tpow3 tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic tasin tacos tcos tatan tsin ttan tsin_cos tagm thypot tfactorial tgamma terf tcbrt tzeta mpf_compat mpfr_compat reuse tsqr tstrtofr tpow_z tget_f tconst_catalan troot tsec tcsc tcot teint tcoth tcsch tsech tstckintc tsubnormal tlngamma tlgamma tzeta_ui tget_ld_2exp tget_set_d64 tj0 tj1 tjn ty0 ty1 tyn tremquo
 
 EXTRA_DIST = tgeneric.c tgeneric_ui.c mpf_compat.h inp_str.data
 
diff --git a/tests/tremquo.c b/tests/tremquo.c
new file mode 100644
index 000000000..2d5814e12
--- /dev/null
+++ b/tests/tremquo.c
@@ -0,0 +1,211 @@
+/* tremquo -- test file for mpfr_remquo and mpfr_remainder
+
+Copyright 2007 Free Software Foundation, Inc.
+Contributed by the Arenaire and Cacao projects, INRIA.
+
+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., 51 Franklin St, Fifth Floor, Boston,
+MA 02110-1301, USA. */
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "mpfr-test.h"
+
+int
+main (int argc, char *argv[])
+{
+  mpfr_t x, y, r;
+  int q[1];
+
+  tests_start_mpfr ();
+
+  mpfr_init (x);
+  mpfr_init (y);
+  mpfr_init (r);
+
+  /* special values */
+  mpfr_set_nan (x);
+  mpfr_set_ui (y, 1, GMP_RNDN);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN(mpfr_nan_p (r));
+
+  mpfr_set_ui (x, 1, GMP_RNDN);
+  mpfr_set_nan (y);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN(mpfr_nan_p (r));
+
+  mpfr_set_inf (x, 1); /* +Inf */
+  mpfr_set_ui (y, 1, GMP_RNDN);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN (mpfr_nan_p (r));
+  
+  mpfr_set_inf (x, 1); /* +Inf */
+  mpfr_set_ui (y, 0, GMP_RNDN);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN (mpfr_nan_p (r));
+  
+  mpfr_set_inf (x, 1); /* +Inf */
+  mpfr_set_inf (y, 1);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN (mpfr_nan_p (r));
+
+  mpfr_set_ui (x, 0, GMP_RNDN);
+  mpfr_set_inf (y, 1);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN (mpfr_cmp_ui (r, 0) == 0 && MPFR_IS_POS (r));
+  MPFR_ASSERTN (q[0] == 0);
+  
+  mpfr_set_ui (x, 0, GMP_RNDN);
+  mpfr_neg (x, x, GMP_RNDN); /* -0 */
+  mpfr_set_inf (y, 1);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN (mpfr_cmp_ui (r, 0) == 0 && MPFR_IS_NEG (r));
+  MPFR_ASSERTN (q[0] == 0);
+  
+  mpfr_set_ui (x, 17, GMP_RNDN);
+  mpfr_set_inf (y, 1);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN (mpfr_cmp (r, x) == 0);
+  MPFR_ASSERTN (q[0] == 0);
+
+  mpfr_set_ui (x, 17, GMP_RNDN);
+  mpfr_set_ui (y, 0, GMP_RNDN);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN (mpfr_nan_p (r));
+  
+  mpfr_set_ui (x, 0, GMP_RNDN);
+  mpfr_set_ui (y, 17, GMP_RNDN);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN (mpfr_cmp_ui (r, 0) == 0 && MPFR_IS_POS (r));
+  MPFR_ASSERTN (q[0] == 0);
+
+  mpfr_set_ui (x, 0, GMP_RNDN);
+  mpfr_neg (x, x, GMP_RNDN);
+  mpfr_set_ui (y, 17, GMP_RNDN);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN (mpfr_cmp_ui (r, 0) == 0 && MPFR_IS_NEG (r));
+  MPFR_ASSERTN (q[0] == 0);
+
+  mpfr_set_prec (x, 53);
+  mpfr_set_prec (y, 53);
+
+  mpfr_set_ui (x, 42, GMP_RNDN);
+  mpfr_set_ui (y, 17, GMP_RNDN);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN (mpfr_cmp_ui (r, 8) == 0);
+  MPFR_ASSERTN (q[0] == 2);
+  
+  mpfr_set_si (x, -42, GMP_RNDN);
+  mpfr_set_ui (y, 17, GMP_RNDN);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN (mpfr_cmp_si (r, -8) == 0);
+  MPFR_ASSERTN (q[0] == -2);
+
+  mpfr_set_prec (x, 100);
+  mpfr_set_prec (y, 50);
+  mpfr_set_ui (x, 42, GMP_RNDN);
+  mpfr_nextabove (x); /* 42 + 2^(-94) */
+  mpfr_set_ui (y, 21, GMP_RNDN);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN (mpfr_cmp_ui_2exp (r, 1, -94) == 0);
+  MPFR_ASSERTN (q[0] == 2);
+  
+  mpfr_set_prec (x, 50);
+  mpfr_set_prec (y, 100);
+  mpfr_set_ui (x, 42, GMP_RNDN);
+  mpfr_nextabove (x); /* 42 + 2^(-44) */
+  mpfr_set_ui (y, 21, GMP_RNDN);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN (mpfr_cmp_ui_2exp (r, 1, -44) == 0);
+  MPFR_ASSERTN (q[0] == 2);
+  
+  mpfr_set_prec (x, 100);
+  mpfr_set_prec (y, 50);
+  mpfr_set_ui (x, 42, GMP_RNDN);
+  mpfr_set_ui (y, 21, GMP_RNDN);
+  mpfr_nextabove (y); /* 21 + 2^(-45) */
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  /* r should be 42 - 2*(21 + 2^(-45)) = -2^(-44) */
+  MPFR_ASSERTN (mpfr_cmp_si_2exp (r, -1, -44) == 0);
+  MPFR_ASSERTN (q[0] == 2);
+
+  mpfr_set_prec (x, 50);
+  mpfr_set_prec (y, 100);
+  mpfr_set_ui (x, 42, GMP_RNDN);
+  mpfr_set_ui (y, 21, GMP_RNDN);
+  mpfr_nextabove (y); /* 21 + 2^(-95) */
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  /* r should be 42 - 2*(21 + 2^(-95)) = -2^(-94) */
+  MPFR_ASSERTN (mpfr_cmp_si_2exp (r, -1, -94) == 0);
+  MPFR_ASSERTN (q[0] == 2);
+
+  /* exercise large quotient */
+  mpfr_set_ui_2exp (x, 1, 65, GMP_RNDN);
+  mpfr_set_ui (y, 1, GMP_RNDN);
+  /* quotient is 2^65, i.e. has 66 bits, and should thus be split on two
+     consecutive limbs */
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  MPFR_ASSERTN (mpfr_cmp_si (r, 0) == 0);
+  MPFR_ASSERTN (q[0] == 0);
+
+  /* another large quotient */
+  mpfr_set_prec (x, 65);
+  mpfr_set_prec (y, 65);
+  mpfr_const_pi (x, GMP_RNDN);
+  mpfr_mul_2exp (x, x, 63, GMP_RNDN);
+  mpfr_const_log2 (y, GMP_RNDN);
+  mpfr_set_prec (r, 10);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  /* q should be 41803643793084085130, r should be 605/2048 */
+  MPFR_ASSERTN (mpfr_cmp_ui_2exp (r, 605, -11) == 0);
+  MPFR_ASSERTN (q[0] == 2);
+
+  /* check cases where quotient is 1.5 +/- eps */
+  mpfr_set_prec (x, 65);
+  mpfr_set_prec (y, 65);
+  mpfr_set_prec (r, 63);
+  mpfr_set_ui (x, 3, GMP_RNDN);
+  mpfr_set_ui (y, 2, GMP_RNDN);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  /* x/y = 1.5, quotient should be 2 (even rule), remainder should be -1 */
+  MPFR_ASSERTN (mpfr_cmp_si (r, -1) == 0);
+  MPFR_ASSERTN (q[0] = 2);
+  mpfr_set_ui (x, 3, GMP_RNDN);
+  mpfr_nextabove (x); /* 3 + 2^(-63) */
+  mpfr_set_ui (y, 2, GMP_RNDN);
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  /* x/y = 1.5 + 2^(-64), quo should be 2, r should be -1 + 2^(-63) */
+  MPFR_ASSERTN (mpfr_add_ui (r, r, 1, GMP_RNDN) == 0);
+  MPFR_ASSERTN (mpfr_cmp_ui_2exp (r, 1, -63) == 0);
+  MPFR_ASSERTN (q[0] = 2);
+  mpfr_set_ui (x, 3, GMP_RNDN);
+  mpfr_set_ui (y, 2, GMP_RNDN);
+  mpfr_nextabove (y); /* 2 + 2^(-63) */
+  mpfr_remquo (q, r, x, y, GMP_RNDN);
+  /* x/y = 1.5 - eps, quo should be 1, r should be 1 - 2^(-63) */
+  MPFR_ASSERTN (mpfr_sub_ui (r, r, 1, GMP_RNDN) == 0);
+  MPFR_ASSERTN (mpfr_cmp_si_2exp (r, -1, -63) == 0);
+  MPFR_ASSERTN (q[0] = 1);
+
+  mpfr_clear (x);
+  mpfr_clear (y);
+  mpfr_clear (r);
+
+  tests_end_mpfr ();
+
+  return 0;
+}