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 } 151