1/* mpfr_subnormalize -- Subnormalize a floating point number 2 emulating sub-normal numbers. 3 4Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 5Contributed by the AriC and Caramel projects, INRIA. 6 7This file is part of the GNU MPFR Library. 8 9The GNU MPFR Library is free software; you can redistribute it and/or modify 10it under the terms of the GNU Lesser General Public License as published by 11the Free Software Foundation; either version 3 of the License, or (at your 12option) any later version. 13 14The GNU MPFR Library is distributed in the hope that it will be useful, but 15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 16or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 17License for more details. 18 19You should have received a copy of the GNU Lesser General Public License 20along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 21http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2251 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 23 24#include "mpfr-impl.h" 25 26/* For MPFR_RNDN, we can have a problem of double rounding. 27 In such a case, this table helps to conclude what to do (y positive): 28 Rounding Bit | Sticky Bit | inexact | Action | new inexact 29 0 | ? | ? | Trunc | sticky 30 1 | 0 | 1 | Trunc | 31 1 | 0 | 0 | Trunc if even | 32 1 | 0 | -1 | AddOneUlp | 33 1 | 1 | ? | AddOneUlp | 34 35 For other rounding mode, there isn't such a problem. 36 Just round it again and merge the ternary values. 37 38 Set the inexact flag if the returned ternary value is non-zero. 39 Set the underflow flag if a second rounding occurred (whether this 40 rounding is exact or not). See 41 https://sympa.inria.fr/sympa/arc/mpfr/2009-06/msg00000.html 42 https://sympa.inria.fr/sympa/arc/mpfr/2009-06/msg00008.html 43 https://sympa.inria.fr/sympa/arc/mpfr/2009-06/msg00010.html 44*/ 45 46int 47mpfr_subnormalize (mpfr_ptr y, int old_inexact, mpfr_rnd_t rnd) 48{ 49 int sign; 50 51 /* The subnormal exponent range is [ emin, emin + MPFR_PREC(y) - 2 ] */ 52 if (MPFR_LIKELY (MPFR_IS_SINGULAR (y) 53 || (MPFR_GET_EXP (y) >= 54 __gmpfr_emin + (mpfr_exp_t) MPFR_PREC (y) - 1))) 55 MPFR_RET (old_inexact); 56 57 mpfr_set_underflow (); 58 sign = MPFR_SIGN (y); 59 60 /* We have to emulate one bit rounding if EXP(y) = emin */ 61 if (MPFR_GET_EXP (y) == __gmpfr_emin) 62 { 63 /* If this is a power of 2, we don't need rounding. 64 It handles cases when |y| = 0.1 * 2^emin */ 65 if (mpfr_powerof2_raw (y)) 66 MPFR_RET (old_inexact); 67 68 /* We keep the same sign for y. 69 Assuming Y is the real value and y the approximation 70 and since y is not a power of 2: 0.5*2^emin < Y < 1*2^emin 71 We also know the direction of the error thanks to ternary value. */ 72 73 if (rnd == MPFR_RNDN) 74 { 75 mp_limb_t *mant, rb ,sb; 76 mp_size_t s; 77 /* We need the rounding bit and the sticky bit. Read them 78 and use the previous table to conclude. */ 79 s = MPFR_LIMB_SIZE (y) - 1; 80 mant = MPFR_MANT (y) + s; 81 rb = *mant & (MPFR_LIMB_HIGHBIT >> 1); 82 if (rb == 0) 83 goto set_min; 84 sb = *mant & ((MPFR_LIMB_HIGHBIT >> 1) - 1); 85 while (sb == 0 && s-- != 0) 86 sb = *--mant; 87 if (sb != 0) 88 goto set_min_p1; 89 /* Rounding bit is 1 and sticky bit is 0. 90 We need to examine old inexact flag to conclude. */ 91 if ((old_inexact > 0 && sign > 0) || 92 (old_inexact < 0 && sign < 0)) 93 goto set_min; 94 /* If inexact != 0, return 0.1*2^(emin+1). 95 Otherwise, rounding bit = 1, sticky bit = 0 and inexact = 0 96 So we have 0.1100000000000000000000000*2^emin exactly. 97 We return 0.1*2^(emin+1) according to the even-rounding 98 rule on subnormals. */ 99 goto set_min_p1; 100 } 101 else if (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG (y))) 102 { 103 set_min: 104 mpfr_setmin (y, __gmpfr_emin); 105 MPFR_RET (-sign); 106 } 107 else 108 { 109 set_min_p1: 110 /* Note: mpfr_setmin will abort if __gmpfr_emax == __gmpfr_emin. */ 111 mpfr_setmin (y, __gmpfr_emin + 1); 112 MPFR_RET (sign); 113 } 114 } 115 else /* Hard case: It is more or less the same problem than mpfr_cache */ 116 { 117 mpfr_t dest; 118 mpfr_prec_t q; 119 int inexact, inex2; 120 121 MPFR_ASSERTD (MPFR_GET_EXP (y) > __gmpfr_emin); 122 123 /* Compute the intermediary precision */ 124 q = (mpfr_uexp_t) MPFR_GET_EXP (y) - __gmpfr_emin + 1; 125 MPFR_ASSERTD (q >= MPFR_PREC_MIN && q < MPFR_PREC (y)); 126 127 /* TODO: perform the rounding in place. */ 128 mpfr_init2 (dest, q); 129 /* Round y in dest */ 130 MPFR_SET_EXP (dest, MPFR_GET_EXP (y)); 131 MPFR_SET_SIGN (dest, sign); 132 MPFR_RNDRAW_EVEN (inexact, dest, 133 MPFR_MANT (y), MPFR_PREC (y), rnd, sign, 134 MPFR_SET_EXP (dest, MPFR_GET_EXP (dest) + 1)); 135 if (MPFR_LIKELY (old_inexact != 0)) 136 { 137 if (MPFR_UNLIKELY (rnd == MPFR_RNDN && 138 (inexact == MPFR_EVEN_INEX || 139 inexact == -MPFR_EVEN_INEX))) 140 { 141 /* if both roundings are in the same direction, we have to go 142 back in the other direction */ 143 if (SAME_SIGN (inexact, old_inexact)) 144 { 145 if (SAME_SIGN (inexact, MPFR_INT_SIGN (y))) 146 mpfr_nexttozero (dest); 147 else 148 mpfr_nexttoinf (dest); 149 inexact = -inexact; 150 } 151 } 152 else if (MPFR_UNLIKELY (inexact == 0)) 153 inexact = old_inexact; 154 } 155 156 inex2 = mpfr_set (y, dest, rnd); 157 MPFR_ASSERTN (inex2 == 0); 158 MPFR_ASSERTN (MPFR_IS_PURE_FP (y)); 159 mpfr_clear (dest); 160 161 MPFR_RET (inexact); 162 } 163} 164