From 1f2c7864507e7a3a9110c0e1d755f7a80df16914 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Wed, 21 Nov 2001 20:07:35 +0000 Subject: mpfr_add_one_ulp and mpfr_sub_one_ulp changed. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@1570 280ebfd0-de03-0410-8827-d642c229c3f4 --- Makefile.am | 2 +- Makefile.in | 38 ++++++++++++++++---------------- add_one_ulp.c | 56 ++++++++++++++++++++++++++++++++++++++++++++++++ add_ulp.c | 69 ----------------------------------------------------------- div_ui.c | 8 +++---- exp.c | 4 ++-- mpfr.h | 4 ++-- mpfr.texi | 23 +++++++++++--------- set_z.c | 69 ++++++++++++++++++++++++++++++++++++++++------------------- sub_one_ulp.c | 63 +++++++++++++++++++++++++++++++++++++++++++++++++++++ tests/tdiv.c | 6 +++--- tests/texp.c | 2 +- tests/tsqrt.c | 2 +- 13 files changed, 213 insertions(+), 133 deletions(-) create mode 100644 add_one_ulp.c delete mode 100644 add_ulp.c create mode 100644 sub_one_ulp.c diff --git a/Makefile.am b/Makefile.am index 24fb080cc..6b7e0fcf0 100644 --- a/Makefile.am +++ b/Makefile.am @@ -5,7 +5,7 @@ include_HEADERS = mpfr.h mpf2mpfr.h mpfr-test.h lib_LIBRARIES = libmpfr.a -libmpfr_a_SOURCES= cputime.h mpfr.h mpfr-impl.h exceptions.c save_expo.c extract.c add.c add1.c add_ui.c add_ulp.c agm.c clear.c cmp.c cmp_abs.c cmp_ui.c div_2exp.c div.c div_ui.c dump.c eq.c exp2.c exp3.c exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c log2.c log.c mul_2exp.c mul.c mul_ui.c neg.c out_str.c pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c rnd_mode.c round.c set.c set_d.c set_dfl_prec.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 trunc.c ui_div.c ui_sub.c urandomb.c sqrtrem.c mpz_set_fr.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 euler.c cos.c sin.c tan.c fma.c hypot.c log1p.c expm1.c log_base_2.c ui_pow.c ui_pow_ui.c log_base_10.c minmax.c dim.c copysign.c gmp_op.c init2.c +libmpfr_a_SOURCES= cputime.h mpfr.h mpfr-impl.h exceptions.c save_expo.c extract.c add.c add1.c add_one_ulp.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_ui.c div_2exp.c div.c div_ui.c dump.c eq.c exp2.c exp3.c exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c log2.c log.c mul_2exp.c mul.c mul_ui.c neg.c out_str.c pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c rnd_mode.c round.c set.c set_d.c set_dfl_prec.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_one_ulp.c sub_ui.c trunc.c ui_div.c ui_sub.c urandomb.c sqrtrem.c mpz_set_fr.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 euler.c cos.c sin.c tan.c fma.c hypot.c log1p.c expm1.c log_base_2.c ui_pow.c ui_pow_ui.c log_base_10.c minmax.c dim.c copysign.c gmp_op.c init2.c libmpfr_a_LIBADD= ceil.o floor.o diff --git a/Makefile.in b/Makefile.in index b7c6117ad..4ee7eb3b9 100644 --- a/Makefile.in +++ b/Makefile.in @@ -75,7 +75,7 @@ include_HEADERS = mpfr.h mpf2mpfr.h mpfr-test.h lib_LIBRARIES = libmpfr.a -libmpfr_a_SOURCES = cputime.h mpfr.h mpfr-impl.h exceptions.c save_expo.c extract.c add.c add1.c add_ui.c add_ulp.c agm.c clear.c cmp.c cmp_abs.c cmp_ui.c div_2exp.c div.c div_ui.c dump.c eq.c exp2.c exp3.c exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c log2.c log.c mul_2exp.c mul.c mul_ui.c neg.c out_str.c pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c rnd_mode.c round.c set.c set_d.c set_dfl_prec.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 trunc.c ui_div.c ui_sub.c urandomb.c sqrtrem.c mpz_set_fr.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 euler.c cos.c sin.c tan.c fma.c hypot.c log1p.c expm1.c log_base_2.c ui_pow.c ui_pow_ui.c log_base_10.c minmax.c dim.c copysign.c gmp_op.c init2.c +libmpfr_a_SOURCES = cputime.h mpfr.h mpfr-impl.h exceptions.c save_expo.c extract.c add.c add1.c add_one_ulp.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_ui.c div_2exp.c div.c div_ui.c dump.c eq.c exp2.c exp3.c exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c log2.c log.c mul_2exp.c mul.c mul_ui.c neg.c out_str.c pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c rnd_mode.c round.c set.c set_d.c set_dfl_prec.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_one_ulp.c sub_ui.c trunc.c ui_div.c ui_sub.c urandomb.c sqrtrem.c mpz_set_fr.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 euler.c cos.c sin.c tan.c fma.c hypot.c log1p.c expm1.c log_base_2.c ui_pow.c ui_pow_ui.c log_base_10.c minmax.c dim.c copysign.c gmp_op.c init2.c libmpfr_a_LIBADD = ceil.o floor.o @@ -93,18 +93,19 @@ CPPFLAGS = @CPPFLAGS@ LIBS = @LIBS@ libmpfr_a_DEPENDENCIES = ceil.o floor.o libmpfr_a_OBJECTS = exceptions.o save_expo.o extract.o add.o add1.o \ -add_ui.o add_ulp.o agm.o clear.o cmp.o cmp_abs.o cmp_ui.o div_2exp.o \ -div.o div_ui.o dump.o eq.o exp2.o exp3.o exp.o get_str.o init.o \ -inp_str.o isinteger.o isinf.o isnan.o isnum.o log2.o log.o mul_2exp.o \ -mul.o mul_ui.o neg.o out_str.o pi.o pow.o pow_si.o pow_ui.o print_raw.o \ -print_rnd_mode.o random2.o random.o reldiff.o rnd_mode.o round.o set.o \ -set_d.o set_dfl_prec.o set_rnd.o set_f.o set_prc_raw.o set_prec.o \ -set_q.o set_si.o set_str.o set_str_raw.o set_ui.o set_z.o sqrt.o \ -sqrt_ui.o sub.o sub1.o sub_ui.o trunc.o ui_div.o ui_sub.o urandomb.o \ -sqrtrem.o mpz_set_fr.o swap.o factorial.o cosh.o sinh.o tanh.o acosh.o \ -asinh.o atanh.o atan.o cmp2.o exp_2.o asin.o euler.o cos.o sin.o tan.o \ -fma.o hypot.o log1p.o expm1.o log_base_2.o ui_pow.o ui_pow_ui.o \ -log_base_10.o minmax.o dim.o copysign.o gmp_op.o init2.o +add_one_ulp.o add_ui.o agm.o clear.o cmp.o cmp_abs.o cmp_ui.o \ +div_2exp.o div.o div_ui.o dump.o eq.o exp2.o exp3.o exp.o get_str.o \ +init.o inp_str.o isinteger.o isinf.o isnan.o isnum.o log2.o log.o \ +mul_2exp.o mul.o mul_ui.o neg.o out_str.o pi.o pow.o pow_si.o pow_ui.o \ +print_raw.o print_rnd_mode.o random2.o random.o reldiff.o rnd_mode.o \ +round.o set.o set_d.o set_dfl_prec.o set_rnd.o set_f.o set_prc_raw.o \ +set_prec.o set_q.o set_si.o set_str.o set_str_raw.o set_ui.o set_z.o \ +sqrt.o sqrt_ui.o sub.o sub1.o sub_one_ulp.o sub_ui.o trunc.o ui_div.o \ +ui_sub.o urandomb.o sqrtrem.o mpz_set_fr.o swap.o factorial.o cosh.o \ +sinh.o tanh.o acosh.o asinh.o atanh.o atan.o cmp2.o exp_2.o asin.o \ +euler.o cos.o sin.o tan.o fma.o hypot.o log1p.o expm1.o log_base_2.o \ +ui_pow.o ui_pow_ui.o log_base_10.o minmax.o dim.o copysign.o gmp_op.o \ +init2.o AR = ar COMPILE = $(CC) $(DEFS) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) CCLD = $(CC) @@ -124,8 +125,8 @@ DISTFILES = $(DIST_COMMON) $(SOURCES) $(HEADERS) $(TEXINFOS) $(EXTRA_DIST) TAR = gtar GZIP_ENV = --best -DEP_FILES = .deps/acosh.P .deps/add.P .deps/add1.P .deps/add_ui.P \ -.deps/add_ulp.P .deps/agm.P .deps/asin.P .deps/asinh.P .deps/atan.P \ +DEP_FILES = .deps/acosh.P .deps/add.P .deps/add1.P .deps/add_one_ulp.P \ +.deps/add_ui.P .deps/agm.P .deps/asin.P .deps/asinh.P .deps/atan.P \ .deps/atanh.P .deps/clear.P .deps/cmp.P .deps/cmp2.P .deps/cmp_abs.P \ .deps/cmp_ui.P .deps/copysign.P .deps/cos.P .deps/cosh.P .deps/dim.P \ .deps/div.P .deps/div_2exp.P .deps/div_ui.P .deps/dump.P .deps/eq.P \ @@ -143,9 +144,10 @@ DEP_FILES = .deps/acosh.P .deps/add.P .deps/add1.P .deps/add_ui.P \ .deps/set_f.P .deps/set_prc_raw.P .deps/set_prec.P .deps/set_q.P \ .deps/set_rnd.P .deps/set_si.P .deps/set_str.P .deps/set_str_raw.P \ .deps/set_ui.P .deps/set_z.P .deps/sin.P .deps/sinh.P .deps/sqrt.P \ -.deps/sqrt_ui.P .deps/sqrtrem.P .deps/sub.P .deps/sub1.P .deps/sub_ui.P \ -.deps/swap.P .deps/tan.P .deps/tanh.P .deps/trunc.P .deps/ui_div.P \ -.deps/ui_pow.P .deps/ui_pow_ui.P .deps/ui_sub.P .deps/urandomb.P +.deps/sqrt_ui.P .deps/sqrtrem.P .deps/sub.P .deps/sub1.P \ +.deps/sub_one_ulp.P .deps/sub_ui.P .deps/swap.P .deps/tan.P \ +.deps/tanh.P .deps/trunc.P .deps/ui_div.P .deps/ui_pow.P \ +.deps/ui_pow_ui.P .deps/ui_sub.P .deps/urandomb.P SOURCES = $(libmpfr_a_SOURCES) OBJECTS = $(libmpfr_a_OBJECTS) diff --git a/add_one_ulp.c b/add_one_ulp.c new file mode 100644 index 000000000..a889ea226 --- /dev/null +++ b/add_one_ulp.c @@ -0,0 +1,56 @@ +/* mpfr_add_one_ulp -- add one unit in last place + +Copyright (C) 1999, 2001 Free Software Foundation, Inc. + +This file is part of the MPFR Library. + +The MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your +option) any later version. + +The MPFR Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "mpfr.h" +#include "mpfr-impl.h" + +/* sets x to x+sign(x)*ulp(x) */ +int +mpfr_add_one_ulp (mpfr_ptr x, mp_rnd_t rnd_mode) +{ + mp_size_t xn; + int sh; + mp_limb_t *xp; + + if (MPFR_IS_NAN(x)) + MPFR_RET_NAN; + + if (MPFR_IS_INF(x) || MPFR_IS_ZERO(x)) + return 0; + + xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB; + sh = xn * BITS_PER_MP_LIMB - MPFR_PREC(x); + xp = MPFR_MANT(x); + if (mpn_add_1 (xp, xp, xn, MP_LIMB_T_ONE << sh)) /* got 1.0000... */ + { + mp_exp_t exp = MPFR_EXP(x); + if (exp == __mpfr_emax) + return mpfr_set_overflow(x, rnd_mode, MPFR_SIGN(x)); + else + { + MPFR_EXP(x)++; + xp[xn-1] = MP_LIMB_T_HIGHBIT; + } + } + return 0; +} diff --git a/add_ulp.c b/add_ulp.c deleted file mode 100644 index a0c0e58f0..000000000 --- a/add_ulp.c +++ /dev/null @@ -1,69 +0,0 @@ -/* mpfr_add_one_ulp, mpfr_sub_one_ulp -- add/subtract one unit in last place - -Copyright (C) 1999, 2001 Free Software Foundation, Inc. - -This file is part of the MPFR Library. - -The MPFR Library is free software; you can redistribute it and/or modify -it under the terms of the GNU Lesser General Public License as published by -the Free Software Foundation; either version 2.1 of the License, or (at your -option) any later version. - -The MPFR Library is distributed in the hope that it will be useful, but -WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public -License for more details. - -You should have received a copy of the GNU Lesser General Public License -along with the MPFR Library; see the file COPYING.LIB. If not, write to -the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, -MA 02111-1307, USA. */ - -#include "gmp.h" -#include "mpfr.h" -#include "mpfr-impl.h" -#include "gmp-impl.h" - -/* sets x to x+sign(x)*2^(MPFR_EXP(x)-MPFR_PREC(x)) */ -void -mpfr_add_one_ulp (mpfr_ptr x) -{ - int xn, sh; - mp_limb_t *xp; - - if (MPFR_IS_INF(x)) - return; - - xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB; - sh = xn * BITS_PER_MP_LIMB - MPFR_PREC(x); - xp = MPFR_MANT(x); - if (mpn_add_1 (xp, xp, xn, MP_LIMB_T_ONE << sh)) /* got 1.0000... */ - { - MPFR_EXP(x)++; - xp[xn-1] = MP_LIMB_T_HIGHBIT; - } - return; -} - -/* sets x to x-sign(x)*ulp(x) */ -void -mpfr_sub_one_ulp(mpfr_ptr x) -{ - int xn, sh; - mp_limb_t *xp; - - if (MPFR_IS_INF(x)) - return; - - xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB; - sh = xn * BITS_PER_MP_LIMB - MPFR_PREC(x); - xp = MPFR_MANT(x); - mpn_sub_1 (xp, xp, xn, MP_LIMB_T_ONE << sh); - if (xp[xn-1] >> (BITS_PER_MP_LIMB - 1) == 0) { - /* was an exact power of two: not normalized any more */ - MPFR_EXP(x)--; - mpn_lshift (xp, xp, xn, 1); - *xp |= MP_LIMB_T_ONE << sh; - } - return; -} diff --git a/div_ui.c b/div_ui.c index c3986f07d..db4aa0c63 100644 --- a/div_ui.c +++ b/div_ui.c @@ -149,12 +149,12 @@ mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) case GMP_RNDU: if (MPFR_SIGN(y) > 0) - mpfr_add_one_ulp (y); + mpfr_add_one_ulp (y, rnd_mode); MPFR_RET(1); /* result is inexact */ case GMP_RNDD: if (MPFR_SIGN(y) < 0) - mpfr_add_one_ulp (y); + mpfr_add_one_ulp (y, rnd_mode); MPFR_RET(-1); /* result is inexact */ case GMP_RNDN: @@ -162,7 +162,7 @@ mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) MPFR_RET(-MPFR_SIGN(x)); else if (sh && d > (MP_LIMB_T_ONE << (sh - 1))) { - mpfr_add_one_ulp (y); + mpfr_add_one_ulp (y, rnd_mode); MPFR_RET(MPFR_SIGN(x)); } else /* sh = 0 or d = 1 << (sh-1) */ @@ -173,7 +173,7 @@ mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) */ if ((sh && inexact) || (!sh && (middle > 0)) || (*yp & (MP_LIMB_T_ONE << sh))) { - mpfr_add_one_ulp (y); + mpfr_add_one_ulp (y, rnd_mode); MPFR_RET(MPFR_SIGN(x)); } else diff --git a/exp.c b/exp.c index c2130880f..089975e6f 100644 --- a/exp.c +++ b/exp.c @@ -107,12 +107,12 @@ mpfr_exp (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_set_ui (y, 1, rnd_mode); if (signx > 0 && rnd_mode == GMP_RNDU) { - mpfr_add_one_ulp (y); + mpfr_add_one_ulp (y, rnd_mode); return 1; } else if (signx < 0 && (rnd_mode == GMP_RNDD || rnd_mode == GMP_RNDZ)) { - mpfr_sub_one_ulp (y); + mpfr_sub_one_ulp (y, rnd_mode); return -1; } return -signx; diff --git a/mpfr.h b/mpfr.h index b887c250c..af7d0a766 100644 --- a/mpfr.h +++ b/mpfr.h @@ -166,7 +166,7 @@ int mpfr_sqrt_ui _PROTO ((mpfr_ptr, unsigned long, mp_rnd_t)); int mpfr_add _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t)); int mpfr_add_ui _PROTO ((mpfr_ptr, mpfr_srcptr, unsigned long, mp_rnd_t)); int mpfr_sub_ui _PROTO ((mpfr_ptr, mpfr_srcptr, unsigned long, mp_rnd_t)); -void mpfr_add_one_ulp _PROTO ((mpfr_ptr)); +int mpfr_add_one_ulp _PROTO ((mpfr_ptr, mp_rnd_t)); int mpfr_sub _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t)); int mpfr_ui_sub _PROTO ((mpfr_ptr, unsigned long, mpfr_srcptr, mp_rnd_t)); void mpfr_reldiff _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t)); @@ -193,7 +193,7 @@ 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)); int mpfr_neg _PROTO((mpfr_ptr, mpfr_srcptr, mp_rnd_t)); -void mpfr_sub_one_ulp _PROTO((mpfr_ptr)); +int mpfr_sub_one_ulp _PROTO((mpfr_ptr, mp_rnd_t)); int mpfr_div_ui _PROTO((mpfr_ptr, mpfr_srcptr, unsigned long int, mp_rnd_t)); int mpfr_ui_div _PROTO((mpfr_ptr, unsigned long int, mpfr_srcptr, mp_rnd_t)); mp_prec_t mpfr_get_prec _PROTO((mpfr_srcptr)); diff --git a/mpfr.texi b/mpfr.texi index 102a20560..389bde992 100644 --- a/mpfr.texi +++ b/mpfr.texi @@ -1308,16 +1308,19 @@ These functions were mainly designed for the implementation of @code{mpfr}, but may be useful for users too. You need to include @code{mpfr-impl.h} to use them. -@deftypefun int mpfr_add_one_ulp (mpfr_t @var{x}) -Add one unit in last place (ulp) to the mantissa of @var{x} -if it is positive or zero, and subtracts one ulp otherwise. -Always return 0 (result is exact). -@end deftypefun - -@deftypefun int mpfr_sub_one_ulp (mpfr_t @var{x}) -Subtract one ulp to @var{x} if it is positive or zero, -and adds one ulp otherwise. -Always return 0 (result is exact). +@deftypefun int mpfr_add_one_ulp (mpfr_t @var{x}, mp_rnd_t @var{rnd}) +Add one unit in last place (ulp) to @var{x} if @var{x} is finite +and positive, subtract one ulp if @var{x} is finite and negative; +otherwise, @var{x} is not changed. +The return value is zero unless an overflow occurs, in which case the +@code{mpfr_add_one_ulp} function behaves like a conventional addition. +@end deftypefun + +@deftypefun int mpfr_sub_one_ulp (mpfr_t @var{x}, mp_rnd_t @var{rnd}) +Subtract one ulp to @var{x} if @var{x} is finite and positive, add one +ulp if @var{x} is finite and negative; otherwise, @var{x} is not changed. +The return value is zero unless an underflow occurs, in which case the +@code{mpfr_sub_one_ulp} function behaves like a conventional subtraction. @end deftypefun @deftypefun int mpfr_can_round (mpfr_t @var{b}, mp_prec_t @var{err}, mp_rnd_t @var{rnd1}, mp_rnd_t @var{rnd2}, mp_prec_t @var{prec}) diff --git a/set_z.c b/set_z.c index 742b806ac..e8a9d55c7 100644 --- a/set_z.c +++ b/set_z.c @@ -20,14 +20,14 @@ the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #include "gmp.h" -#include "mpfr.h" -#include "mpfr-impl.h" #include "gmp-impl.h" #include "longlong.h" +#include "mpfr.h" +#include "mpfr-impl.h" /* set f to the integer z */ int -mpfr_set_z (mpfr_ptr f, mpz_srcptr z, mp_rnd_t rnd) +mpfr_set_z (mpfr_ptr f, mpz_srcptr z, mp_rnd_t rnd_mode) { mp_size_t fn, zn, dif, sh; int k, sign_z; @@ -58,7 +58,7 @@ mpfr_set_z (mpfr_ptr f, mpz_srcptr z, mp_rnd_t rnd) sh = fn*BITS_PER_MP_LIMB-MPFR_PREC(f); cc = *fp & ((MP_LIMB_T_ONE << sh) - 1); *fp = *fp & ~cc; - if (rnd==GMP_RNDN) { + if (rnd_mode == GMP_RNDN) { if (sh) c2 = MP_LIMB_T_ONE << (sh - 1); else { /* sh=0 */ c2 = MP_LIMB_T_HIGHBIT; @@ -67,25 +67,50 @@ mpfr_set_z (mpfr_ptr f, mpz_srcptr z, mp_rnd_t rnd) if (dif>0 && k) cc += zp[dif-1] >> (BITS_PER_MP_LIMB-k); } /* now compares cc to c2 */ - if (cc>c2) { mpfr_add_one_ulp(f); return cc; } - else if (cc0 && (cc=zp[dif-1])==0) dif--; - if (cc) { mpfr_add_one_ulp(f); return cc; } - else /* exactly in middle: inexact in both cases */ - if (*fp & (MP_LIMB_T_ONE << sh)) { mpfr_add_one_ulp(f); return 1; } - else return 1; - } - } - else if ((sign_z>0 && rnd==GMP_RNDU) || (sign_z<0 && rnd==GMP_RNDD)) { - /* round towards infinity */ - /* result is exact iff all remaining bits are zero */ - if (dif>0 && cc==0) cc=zp[--dif]<0) cc=zp[--dif]; - if (cc) { mpfr_add_one_ulp(f); return 1; } - else return 0; + if (cc > c2) + { + mpfr_add_one_ulp(f, rnd_mode); + return cc; + } + else if (cc < c2) + goto towards_zero; + else + { + cc = 0; + while (dif > 0 && (cc = zp[dif-1]) == 0) + dif--; + if (cc) + { + mpfr_add_one_ulp(f, rnd_mode); + return cc; + } + else /* exactly in middle: inexact in both cases */ + if (*fp & (MP_LIMB_T_ONE << sh)) + { + mpfr_add_one_ulp(f, rnd_mode); + return 1; + } + else + return 1; + } } + else if ((sign_z > 0 && rnd_mode == GMP_RNDU) || + (sign_z < 0 && rnd_mode == GMP_RNDD)) + { + /* round towards infinity */ + /* result is exact iff all remaining bits are zero */ + if (dif > 0 && cc == 0) + cc = zp[--dif] << k; + while (cc == 0 && dif > 0) + cc = zp[--dif]; + if (cc) + { + mpfr_add_one_ulp(f, rnd_mode); + return 1; + } + else + return 0; + } else { /* round towards zero */ /* result is exact iff all remaining bits are zero */ towards_zero: diff --git a/sub_one_ulp.c b/sub_one_ulp.c new file mode 100644 index 000000000..d4ed213a1 --- /dev/null +++ b/sub_one_ulp.c @@ -0,0 +1,63 @@ +/* mpfr_sub_one_ulp -- subtract one unit in last place + +Copyright (C) 1999, 2001 Free Software Foundation, Inc. + +This file is part of the MPFR Library. + +The MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your +option) any later version. + +The MPFR Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "mpfr.h" +#include "mpfr-impl.h" + +/* sets x to x-sign(x)*ulp(x) */ +int +mpfr_sub_one_ulp(mpfr_ptr x, mp_rnd_t rnd_mode) +{ + mp_size_t xn; + int sh; + mp_limb_t *xp; + + if (MPFR_IS_NAN(x)) + MPFR_RET_NAN; + + if (MPFR_IS_INF(x) || MPFR_IS_ZERO(x)) + return 0; + + if (MPFR_PREC(x) == 1) + return mpfr_sub(x, x, x, rnd_mode); /* to generate the correct zero */ + + xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB; + sh = xn * BITS_PER_MP_LIMB - MPFR_PREC(x); + xp = MPFR_MANT(x); + mpn_sub_1 (xp, xp, xn, MP_LIMB_T_ONE << sh); + if (xp[xn-1] >> (BITS_PER_MP_LIMB - 1) == 0) + { /* was an exact power of two: not normalized any more */ + mp_exp_t exp = MPFR_EXP(x); + if (exp == __mpfr_emin) + return mpfr_set_underflow(x, rnd_mode, MPFR_SIGN(x)); + else + { + int i; + MPFR_EXP(x)--; + xp[0] = (sh + 1 == BITS_PER_MP_LIMB) ? 0 : MP_LIMB_T_MAX << (sh + 1); + for (i = 1; i < xn; i++) + xp[i] = MP_LIMB_T_MAX; + } + } + return 0; +} diff --git a/tests/tdiv.c b/tests/tdiv.c index 9ce84c61c..94a5a4de8 100644 --- a/tests/tdiv.c +++ b/tests/tdiv.c @@ -216,7 +216,7 @@ void check_lowr () } else if (c == 2) { - mpfr_add_one_ulp(z); + mpfr_add_one_ulp(z, GMP_RNDN); if (mpfr_cmp(z2, z)) { fprintf(stderr, "Error in mpfr_div [even rnd?] rnd=GMP_RNDN\n"); @@ -229,7 +229,7 @@ void check_lowr () } else if (c == -2) { - mpfr_sub_one_ulp(z); + mpfr_sub_one_ulp(z, GMP_RNDN); if (mpfr_cmp(z2, z)) { fprintf(stderr, "Error in mpfr_div [even rnd?] rnd=GMP_RNDN\n"); @@ -256,7 +256,7 @@ void check_lowr () mpfr_random(tmp); mpfr_mul(x, z, tmp, GMP_RNDN); mpfr_set(y, tmp, GMP_RNDD); - mpfr_add_one_ulp(x); + mpfr_add_one_ulp(x, GMP_RNDN); c = mpfr_div(z2, x, y, GMP_RNDD); mpfr_div(z3, x, y, GMP_RNDD); diff --git a/tests/texp.c b/tests/texp.c index 129b394ef..5d5362db0 100644 --- a/tests/texp.c +++ b/tests/texp.c @@ -121,7 +121,7 @@ int check_worst_case (double X, double expx) mpfr_exp(x, x, GMP_RNDN); mpfr_set_d(x, X, GMP_RNDN); mpfr_exp(x, x, GMP_RNDU); - mpfr_add_one_ulp(y); + mpfr_add_one_ulp(y, GMP_RNDN); if (mpfr_cmp(x,y)) { fprintf(stderr, "exp(x) rounded towards +infinity is wrong\n"); exit(1); } diff --git a/tests/tsqrt.c b/tests/tsqrt.c index 25cc6d54e..9be24323c 100644 --- a/tests/tsqrt.c +++ b/tests/tsqrt.c @@ -188,7 +188,7 @@ void special () { mpfr_set_prec (z, p); mpfr_set_ui (z, 1, GMP_RNDN); - mpfr_add_one_ulp (z); + mpfr_add_one_ulp (z, GMP_RNDN); mpfr_sqrt (x, z, GMP_RNDU); if (mpfr_cmp_ui (x, 2)) { -- cgit v1.2.1