From f30e1255c22e3d1beca6249e5d5946b574d8e0b1 Mon Sep 17 00:00:00 2001 From: zimmerma Date: Wed, 2 May 2007 10:12:11 +0000 Subject: added mpfr_remquo and mpfr_remainder git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@4420 280ebfd0-de03-0410-8827-d642c229c3f4 --- Makefile.am | 2 +- TODO | 13 ---- mpfr.h | 4 ++ mpfr.texi | 13 ++++ remquo.c | 208 +++++++++++++++++++++++++++++++++++++++++++++++++++++ tests/Makefile.am | 2 +- tests/tremquo.c | 211 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 7 files changed, 438 insertions(+), 15 deletions(-) create mode 100644 remquo.c create mode 100644 tests/tremquo.c 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 , 29 Nov 2003] - 1/sqrt(x) [Regis Dupont , 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 +#include + +#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; +} -- cgit v1.2.1