As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option.
Alternatively you can here view or download the uninterpreted source code file.
For more information about "subnormal.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 4.0.1_vs_4.0.2.

1 /* mpfr_subnormalize -- Subnormalize a floating point number 2 emulating sub-normal numbers. 3 4 Copyright 2005-2019 Free Software Foundation, Inc. 5 Contributed by the AriC and Caramba 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 https://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 /* 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 46 int 47 mpfr_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 as 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 /* subnormal range, thus no overflow */ 148 { 149 mpfr_nexttoinf (dest); 150 MPFR_ASSERTD(!MPFR_IS_INF (dest)); 151 } 152 inexact = -inexact; 153 } 154 } 155 else if (MPFR_UNLIKELY (inexact == 0)) 156 inexact = old_inexact; 157 } 158 159 inex2 = mpfr_set (y, dest, rnd); 160 MPFR_ASSERTN (inex2 == 0); 161 MPFR_ASSERTN (MPFR_IS_PURE_FP (y)); 162 mpfr_clear (dest); 163 164 MPFR_RET (inexact); 165 } 166 }