contrib/mpfr/src/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, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
   5 Contributed by the AriC and Caramel 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 3 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.LESSER.  If not, see
  21 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  22 51 Franklin St, Fifth Floor, Boston, 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_setmax (x, __gmpfr_emax);
  32       return;
  33     }
  34   else if (MPFR_UNLIKELY( MPFR_IS_ZERO(x) ))
  35     {
  36       MPFR_CHANGE_SIGN(x);
  37       mpfr_setmin (x, __gmpfr_emin);
  38     }
  39   else
  40     {
  41       mp_size_t xn;
  42       int sh;
  43       mp_limb_t *xp;
  44
  45       xn = MPFR_LIMB_SIZE (x);
  46       MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC(x));
  47       xp = MPFR_MANT(x);
  48       mpn_sub_1 (xp, xp, xn, MPFR_LIMB_ONE << sh);
  49       if (MPFR_UNLIKELY( MPFR_LIMB_MSB(xp[xn-1]) == 0) )
  50         { /* was an exact power of two: not normalized any more */
  51           mpfr_exp_t exp = MPFR_EXP (x);
  52           if (MPFR_UNLIKELY(exp == __gmpfr_emin))
  53             MPFR_SET_ZERO(x);
  54           else
  55             {
  56               mp_size_t i;
  57               MPFR_SET_EXP (x, exp - 1);
  58               xp[0] = MP_LIMB_T_MAX << sh;
  59               for (i = 1; i < xn; i++)
  60                 xp[i] = MP_LIMB_T_MAX;
  61             }
  62         }
  63     }
  64 }
  65
  66 void
  67 mpfr_nexttoinf (mpfr_ptr x)
  68 {
  69   if (MPFR_UNLIKELY(MPFR_IS_INF(x)))
  70     return;
  71   else if (MPFR_UNLIKELY(MPFR_IS_ZERO(x)))
  72     mpfr_setmin (x, __gmpfr_emin);
  73   else
  74     {
  75       mp_size_t xn;
  76       int sh;
  77       mp_limb_t *xp;
  78
  79       xn = MPFR_LIMB_SIZE (x);
  80       MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC(x));
  81       xp = MPFR_MANT(x);
  82       if (MPFR_UNLIKELY( mpn_add_1 (xp, xp, xn, MPFR_LIMB_ONE << sh)) )
  83         /* got 1.0000... */
  84         {
  85           mpfr_exp_t exp = MPFR_EXP (x);
  86           if (MPFR_UNLIKELY(exp == __gmpfr_emax))
  87             MPFR_SET_INF(x);
  88           else
  89             {
  90               MPFR_SET_EXP (x, exp + 1);
  91               xp[xn-1] = MPFR_LIMB_HIGHBIT;
  92             }
  93         }
  94     }
  95 }
  96
  97 void
  98 mpfr_nextabove (mpfr_ptr x)
  99 {
 100   if (MPFR_UNLIKELY(MPFR_IS_NAN(x)))
 101     {
 102       __gmpfr_flags |= MPFR_FLAGS_NAN;
 103       return;
 104     }
 105   if (MPFR_IS_NEG(x))
 106     mpfr_nexttozero (x);
 107   else
 108     mpfr_nexttoinf (x);
 109 }
 110
 111 void
 112 mpfr_nextbelow (mpfr_ptr x)
 113 {
 114   if (MPFR_UNLIKELY(MPFR_IS_NAN(x)))
 115     {
 116       __gmpfr_flags |= MPFR_FLAGS_NAN;
 117       return;
 118     }
 119
 120   if (MPFR_IS_NEG(x))
 121     mpfr_nexttoinf (x);
 122   else
 123     mpfr_nexttozero (x);
 124 }
 125
 126 void
 127 mpfr_nexttoward (mpfr_ptr x, mpfr_srcptr y)
 128 {
 129   int s;
 130
 131   if (MPFR_UNLIKELY(MPFR_IS_NAN(x)))
 132     {
 133       __gmpfr_flags |= MPFR_FLAGS_NAN;
 134       return;
 135     }
 136   else if (MPFR_UNLIKELY(MPFR_IS_NAN(x) || MPFR_IS_NAN(y)))
 137     {
 138       MPFR_SET_NAN(x);
 139       __gmpfr_flags |= MPFR_FLAGS_NAN;
 140       return;
 141     }
 142
 143   s = mpfr_cmp (x, y);
 144   if (s == 0)
 145     return;
 146   else if (s < 0)
 147     mpfr_nextabove (x);
 148   else
 149     mpfr_nextbelow (x);
 150 }