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

12 files changed, 179 insertions, 99 deletions

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_ulp.c b/add_one_ulp.c
index a0c0e58f0..a889ea226 100644
--- a/add_ulp.c
+++ b/add_one_ulp.c
@@ -1,4 +1,4 @@
-/* mpfr_add_one_ulp, mpfr_sub_one_ulp -- add/subtract one unit in last place
+/* mpfr_add_one_ulp -- add one unit in last place
 
 Copyright (C) 1999, 2001 Free Software Foundation, Inc.
 
@@ -20,50 +20,37 @@ 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"
-#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)
+/* sets x to x+sign(x)*ulp(x) */
+int
+mpfr_add_one_ulp (mpfr_ptr x, mp_rnd_t rnd_mode)
 {
-  int xn, sh;
+  mp_size_t xn;
+  int sh;
   mp_limb_t *xp;
 
-  if (MPFR_IS_INF(x))
-    return;
+  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... */
     {
-      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;
+      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/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 (cc<c2) goto towards_zero;
-      else {
-	cc=0;
-	while (dif>0 && (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]<<k;
-      while (cc==0 && 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))
 	{