contrib/mpfr/next.c

   1 /* mpfr_nextabove, mpfr_nextbelow, mpfr_nexttoward -- next representable
   2 floating-point number
   3
   4 Copyright 1999, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   5 Contributed by the Arenaire and Cacao projects, INRIA.
   6
   7 This file is part of the GNU MPFR Library.
   8
   9 The GNU MPFR Library is free software; you can redistribute it and/or modify
  10 it under the terms of the GNU Lesser General Public License as published by
  11 the Free Software Foundation; either version 2.1 of the License, or (at your
  12 option) any later version.
  13
  14 The GNU MPFR Library is distributed in the hope that it will be useful, but
  15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  16 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  17 License for more details.
  18
  19 You should have received a copy of the GNU Lesser General Public License
  20 along with the GNU MPFR Library; see the file COPYING.LIB.  If not, write to
  21 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  22 MA 02110-1301, USA. */
  23
  24 #include "mpfr-impl.h"
  25
  26 void
  27 mpfr_nexttozero (mpfr_ptr x)
  28 {
  29   if (MPFR_UNLIKELY(MPFR_IS_INF(x)))
  30     {
  31       MPFR_CLEAR_FLAGS(x);
  32       mpfr_setmax (x, __gmpfr_emax);
  33       return;
  34     }
  35   else if (MPFR_UNLIKELY( MPFR_IS_ZERO(x) ))
  36     {
  37       MPFR_CHANGE_SIGN(x);
  38       mpfr_setmin (x, __gmpfr_emin);
  39     }
  40   else
  41     {
  42       mp_size_t xn;
  43       int sh;
  44       mp_limb_t *xp;
  45
  46       xn = MPFR_LIMB_SIZE (x);
  47       MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC(x));
  48       xp = MPFR_MANT(x);
  49       mpn_sub_1 (xp, xp, xn, MPFR_LIMB_ONE << sh);
  50       if (MPFR_UNLIKELY( MPFR_LIMB_MSB(xp[xn-1]) == 0) )
  51         { /* was an exact power of two: not normalized any more */
  52           mp_exp_t exp = MPFR_EXP (x);
  53           if (MPFR_UNLIKELY(exp == __gmpfr_emin))
  54             MPFR_SET_ZERO(x);
  55           else
  56             {
  57               mp_size_t i;
  58               MPFR_SET_EXP (x, exp - 1);
  59               xp[0] = MP_LIMB_T_MAX << sh;
  60               for (i = 1; i < xn; i++)
  61                 xp[i] = MP_LIMB_T_MAX;
  62             }
  63         }
  64     }
  65 }
  66
  67 void
  68 mpfr_nexttoinf (mpfr_ptr x)
  69 {
  70   if (MPFR_UNLIKELY(MPFR_IS_INF(x)))
  71     return;
  72   else if (MPFR_UNLIKELY(MPFR_IS_ZERO(x)))
  73     mpfr_setmin (x, __gmpfr_emin);
  74   else
  75     {
  76       mp_size_t xn;
  77       int sh;
  78       mp_limb_t *xp;
  79
  80       xn = MPFR_LIMB_SIZE (x);
  81       MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC(x));
  82       xp = MPFR_MANT(x);
  83       if (MPFR_UNLIKELY( mpn_add_1 (xp, xp, xn, MPFR_LIMB_ONE << sh)) )
  84         /* got 1.0000... */
  85         {
  86           mp_exp_t exp = MPFR_EXP (x);
  87           if (MPFR_UNLIKELY(exp == __gmpfr_emax))
  88             MPFR_SET_INF(x);
  89           else
  90             {
  91               MPFR_SET_EXP (x, exp + 1);
  92               xp[xn-1] = MPFR_LIMB_HIGHBIT;
  93             }
  94         }
  95     }
  96 }
  97
  98 void
  99 mpfr_nextabove (mpfr_ptr x)
 100 {
 101   if (MPFR_UNLIKELY(MPFR_IS_NAN(x)))
 102     {
 103       __gmpfr_flags |= MPFR_FLAGS_NAN;
 104       return;
 105     }
 106   if (MPFR_IS_NEG(x))
 107     mpfr_nexttozero (x);
 108   else
 109     mpfr_nexttoinf (x);
 110 }
 111
 112 void
 113 mpfr_nextbelow (mpfr_ptr x)
 114 {
 115   if (MPFR_UNLIKELY(MPFR_IS_NAN(x)))
 116     {
 117       __gmpfr_flags |= MPFR_FLAGS_NAN;
 118       return;
 119     }
 120
 121   if (MPFR_IS_NEG(x))
 122     mpfr_nexttoinf (x);
 123   else
 124     mpfr_nexttozero (x);
 125 }
 126
 127 void
 128 mpfr_nexttoward (mpfr_ptr x, mpfr_srcptr y)
 129 {
 130   int s;
 131
 132   if (MPFR_UNLIKELY(MPFR_IS_NAN(x)))
 133     {
 134       __gmpfr_flags |= MPFR_FLAGS_NAN;
 135       return;
 136     }
 137   else if (MPFR_UNLIKELY(MPFR_IS_NAN(x) || MPFR_IS_NAN(y)))
 138     {
 139       MPFR_SET_NAN(x);
 140       __gmpfr_flags |= MPFR_FLAGS_NAN;
 141       return;
 142     }
 143
 144   s = mpfr_cmp (x, y);
 145   if (s == 0)
 146     return;
 147   else if (s < 0)
 148     mpfr_nextabove (x);
 149   else
 150     mpfr_nextbelow (x);
 151 }