New internal functions mpfr_setmin and mpfr_setmax.

New functions mpfr_nextabove, mpfr_nextbelow, mpfr_nexttoward.
Small fix in mpfr_sub_one_ulp.


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

10 files changed, 253 insertions, 23 deletions

diff --git a/Makefile.am b/Makefile.am
index fc38679d6..da30ea598 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -5,7 +5,7 @@ include_HEADERS = mpfr.h mpf2mpfr.h
 
 lib_LIBRARIES = libmpfr.a
 
-libmpfr_a_SOURCES = mpfr.h mpf2mpfr.h mpfr-impl.h mpfr-math.h mpfr-test.h exceptions.c save_expo.c extract.c uceil_exp2.c uceil_log2.c ufloor_log2.c add.c add1.c add_one_ulp.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_si.c cmp_ui.c div_2exp.c div_2si.c div_2ui.c div.c div_ui.c dump.c eq.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 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 rnd_mode.c round_prec.c set.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_one_ulp.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
+libmpfr_a_SOURCES = mpfr.h mpf2mpfr.h mpfr-impl.h mpfr-math.h mpfr-test.h exceptions.c save_expo.c extract.c uceil_exp2.c uceil_log2.c ufloor_log2.c add.c add1.c add_one_ulp.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_si.c cmp_ui.c div_2exp.c div_2si.c div_2ui.c div.c div_ui.c dump.c eq.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 rnd_mode.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_one_ulp.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
 
 libmpfr_a_LIBADD = @LIBOBJS@
 
diff --git a/TODO b/TODO
index be801c6fe..341f32edc 100644
--- a/TODO
+++ b/TODO
@@ -35,8 +35,6 @@ New functions to implement:
 - mpz_set_fr (to set a mpz from a mpfr, with a rounding mode)
 - those from LIA: missing secant, cosecant, cotangent (trigo/hyperbolic)
 - atan2
-- nextafter/nextforward : X + epsilon if X < Y, X - epsilon if X > Y
-- nextabove/nextbelow ?
 - mpfr_greater_p (Kevin, Vincent ok)
                          =     <     >     unordered
 mpfr_greater_p           0     0     1     0
@@ -102,3 +100,5 @@ Miscellaneous:
 - make MPFR thread-safe (e.g. temporary changes of emin/emax in mpfr_set_q)
 
 - randomized tests: do like GMP (see misc.c), but using MPFR_CHECK_RANDOMIZE
+
+- add tests for nextabove/nextbelow/nexttoward.
diff --git a/exceptions.c b/exceptions.c
index ae628551f..be922a08c 100644
--- a/exceptions.c
+++ b/exceptions.c
@@ -196,14 +196,7 @@ mpfr_set_underflow (mpfr_ptr x, mp_rnd_t rnd_mode, int sign)
       || (rnd_mode == GMP_RNDU && sign > 0)
       || (rnd_mode == GMP_RNDD && sign < 0))
     {
-      mp_size_t xn;
-      mp_limb_t *xp;
-
-      MPFR_EXP(x) = __mpfr_emin;
-      xn = (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB;
-      xp = MPFR_MANT(x);
-      xp[xn] = MPFR_LIMB_HIGHBIT;
-      MPN_ZERO(xp, xn);
+      mpfr_setmin (x);
       inex = 1;
     }
   else
@@ -228,17 +221,7 @@ mpfr_set_overflow (mpfr_ptr x, mp_rnd_t rnd_mode, int sign)
   if ((rnd_mode == GMP_RNDU && sign < 0)
    || (rnd_mode == GMP_RNDD && sign > 0))
     {
-      mp_size_t xn, i;
-      int sh;
-      mp_limb_t *xp;
-
-      MPFR_EXP(x) = __mpfr_emax;
-      xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB;
-      sh = (mp_prec_t) xn * BITS_PER_MP_LIMB - MPFR_PREC(x);
-      xp = MPFR_MANT(x);
-      xp[0] = MP_LIMB_T_MAX << sh;
-      for (i = 1; i < xn; i++)
-        xp[i] = MP_LIMB_T_MAX;
+      mpfr_setmax (x);
       inex = -1;
     }
   else
diff --git a/mpfr-impl.h b/mpfr-impl.h
index 1d861d4dc..d0b2fad1d 100644
--- a/mpfr-impl.h
+++ b/mpfr-impl.h
@@ -198,6 +198,8 @@ unsigned long _mpfr_cuberoot _PROTO ((unsigned long));
 int mpfr_exp_2 _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
 int mpfr_exp3 _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
 int mpfr_powerof2_raw _PROTO ((mpfr_srcptr));
+void mpfr_setmax _PROTO ((mpfr_ptr));
+void mpfr_setmin _PROTO ((mpfr_ptr));
 
 #define mpfr_round_raw(yp, xp, xprec, neg, yprec, r, inexp) \
   mpfr_round_raw_generic((yp), (xp), (xprec), (neg), (yprec), (r), (inexp), 0)
diff --git a/mpfr.h b/mpfr.h
index 24cdef9d7..7fd824b65 100644
--- a/mpfr.h
+++ b/mpfr.h
@@ -159,6 +159,9 @@ void mpfr_random _PROTO ((mpfr_ptr));
 void mpfr_random2 _PROTO ((mpfr_ptr, mp_size_t, mp_exp_t)); 
 void mpfr_urandomb _PROTO ((mpfr_ptr, gmp_randstate_t)); 
 void mpfr_clear _PROTO ((mpfr_ptr)); 
+void mpfr_nextabove _PROTO ((mpfr_ptr));
+void mpfr_nextbelow _PROTO ((mpfr_ptr));
+void mpfr_nexttoward _PROTO ((mpfr_ptr, mpfr_srcptr));
 void mpfr_set_str_raw _PROTO ((mpfr_ptr, char *));
 int mpfr_set_str _PROTO ((mpfr_ptr, __gmp_const char *, int, mp_rnd_t));
 int mpfr_init_set_str _PROTO ((mpfr_ptr, char *, int, mp_rnd_t));
diff --git a/mpfr.texi b/mpfr.texi
index ad689a41c..da9946531 100644
--- a/mpfr.texi
+++ b/mpfr.texi
@@ -1347,6 +1347,24 @@ The returned value is zero when the result is exact, positive when it is
 greater than the exact value, and negative when it is smaller.
 @end deftypefun
 
+@deftypefun void mpfr_nexttoward (mpfr_t @var{x}, mpfr_t @var{y})
+If @var{x} or @var{y} is NaN, set @var{x} to NaN; the global invalid
+flag is set too. Otherwise, if @var{x} is different from @var{y},
+replace @var{x} by the next floating-point number (with the precision
+of @var{x} and the current exponent range) in the direction of @var{y}
+(the infinite values are seen as the smallest and largest floating-point
+numbers). If the result is zero, it keeps the same sign. No underflow or
+overflow is generated.
+@end deftypefun
+
+@deftypefun void mpfr_nextabove (mpfr_t @var{x})
+Equivalent to @code{mpfr_nexttoward} where @var{y} is plus infinity.
+@end deftypefun
+
+@deftypefun void mpfr_nextbelow (mpfr_t @var{x})
+Equivalent to @code{mpfr_nexttoward} where @var{y} is minus infinity.
+@end deftypefun
+
 @deftypefun void mpfr_urandomb (mpfr_t @var{rop}, gmp_randstate_t @var{state})
 Generate a uniformly distributed random float in the interval 0 <= X < 1.
 @end deftypefun
diff --git a/next.c b/next.c
new file mode 100644
index 000000000..f545ee964
--- /dev/null
+++ b/next.c
@@ -0,0 +1,141 @@
+/* mpfr_nextabove, mpfr_nextbelow, mpfr_nexttoward -- next representable
+floating-point number
+
+Copyright 1999, 2001, 2002 Free Software Foundation.
+Contributed by the Spaces project, INRIA Lorraine.
+
+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"
+
+static void
+mpfr_nexttozero (mpfr_ptr x)
+{
+  if (MPFR_IS_INF(x))
+    {
+      MPFR_CLEAR_FLAGS(x);
+      mpfr_setmax (x);
+      return;
+    }
+
+  if (MPFR_IS_ZERO(x))
+    {
+      MPFR_CHANGE_SIGN(x);
+      mpfr_setmin (x);
+    }
+  else
+    {
+      mp_size_t xn;
+      int sh;
+      mp_limb_t *xp;
+
+      xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB;
+      sh = (mp_prec_t) 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)
+            MPFR_SET_ZERO(x);
+          else
+            {
+              mp_size_t i;
+              MPFR_EXP(x)--;
+              xp[0] = MP_LIMB_T_MAX << sh;
+              for (i = 1; i < xn; i++)
+                xp[i] = MP_LIMB_T_MAX;
+            }
+        }
+    }
+}
+
+static void
+mpfr_nexttoinf (mpfr_ptr x)
+{
+  if (MPFR_IS_INF(x))
+    return;
+
+  if (MPFR_IS_ZERO(x))
+    mpfr_setmin (x);
+  else
+    {
+      mp_size_t xn;
+      int sh;
+      mp_limb_t *xp;
+
+      xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB;
+      sh = (mp_prec_t) 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)
+            MPFR_SET_INF(x);
+          else
+            {
+              MPFR_EXP(x)++;
+              xp[xn-1] = MPFR_LIMB_HIGHBIT;
+            }
+        }
+    }
+}
+
+void
+mpfr_nextabove (mpfr_ptr x)
+{
+  if (MPFR_IS_NAN(x))
+    MPFR_RET_NAN;
+
+  if (MPFR_SIGN(x) < 0)
+    mpfr_nexttozero (x);
+  else
+    mpfr_nexttoinf (x);
+}
+
+void
+mpfr_nextbelow (mpfr_ptr x)
+{
+  if (MPFR_IS_NAN(x))
+    MPFR_RET_NAN;
+
+  if (MPFR_SIGN(x) < 0)
+    mpfr_nexttoinf (x);
+  else
+    mpfr_nexttozero (x);
+}
+
+void
+mpfr_nexttoward (mpfr_ptr x, mpfr_srcptr y)
+{
+  int s;
+
+  if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y))
+    MPFR_RET_NAN;
+
+  s = mpfr_cmp (x, y);
+  if (s == 0)
+    return;
+  if (s < 0)
+    mpfr_nextabove (x);
+  else
+    mpfr_nextbelow (x);
+}
diff --git a/setmax.c b/setmax.c
new file mode 100644
index 000000000..4e9bf8d8e
--- /dev/null
+++ b/setmax.c
@@ -0,0 +1,43 @@
+/* mpfr_setmax -- maximum representable floating-point number (raw version)
+
+Copyright 2002 Free Software Foundation.
+Contributed by the Spaces project, INRIA Lorraine.
+
+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"
+
+/* Note: the flags are not cleared and the current sign is kept. */
+
+void mpfr_setmax (mpfr_ptr x)
+{
+  mp_size_t xn, i;
+  int sh;
+  mp_limb_t *xp;
+
+  MPFR_EXP(x) = __mpfr_emax;
+  xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB;
+  sh = (mp_prec_t) xn * BITS_PER_MP_LIMB - MPFR_PREC(x);
+  xp = MPFR_MANT(x);
+  xp[0] = MP_LIMB_T_MAX << sh;
+  for (i = 1; i < xn; i++)
+    xp[i] = MP_LIMB_T_MAX;
+}
diff --git a/setmin.c b/setmin.c
new file mode 100644
index 000000000..42eb30693
--- /dev/null
+++ b/setmin.c
@@ -0,0 +1,40 @@
+/* mpfr_setmin -- minimum representable floating-point number (raw version)
+
+Copyright 2002 Free Software Foundation.
+Contributed by the Spaces project, INRIA Lorraine.
+
+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"
+
+/* Note: the flags are not cleared and the current sign is kept. */
+
+void mpfr_setmin (mpfr_ptr x)
+{
+  mp_size_t xn;
+  mp_limb_t *xp;
+
+  MPFR_EXP(x) = __mpfr_emin;
+  xn = (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB;
+  xp = MPFR_MANT(x);
+  xp[xn] = MPFR_LIMB_HIGHBIT;
+  MPN_ZERO(xp, xn);
+}
diff --git a/sub_one_ulp.c b/sub_one_ulp.c
index 95dcbd073..323ab4fd1 100644
--- a/sub_one_ulp.c
+++ b/sub_one_ulp.c
@@ -53,7 +53,7 @@ mpfr_sub_one_ulp(mpfr_ptr x, mp_rnd_t rnd_mode)
         return mpfr_set_underflow(x, rnd_mode, MPFR_SIGN(x));
       else
         {
-          int i;
+          mp_size_t 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++)