/* mpfr_subnormalize -- Subnormalize a floating point number emulating sub-normal numbers. Copyright 2005-2021 Free Software Foundation, Inc. Contributed by the AriC and Caramba projects, INRIA. This file is part of the GNU MPFR Library. The GNU MPFR Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. The GNU MPFR Library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ #include "mpfr-impl.h" /* For MPFR_RNDN, we can have a problem of double rounding. In such a case, this table helps to conclude what to do (y positive): Rounding Bit | Sticky Bit | inexact | Action | new inexact 0 | ? | ? | Trunc | sticky 1 | 0 | 1 | Trunc | 1 | 0 | 0 | Trunc if even | 1 | 0 | -1 | AddOneUlp | 1 | 1 | ? | AddOneUlp | For other rounding mode, there isn't such a problem. Just round it again and merge the ternary values. Set the inexact flag if the returned ternary value is non-zero. Set the underflow flag if a second rounding occurred (whether this rounding is exact or not). See https://sympa.inria.fr/sympa/arc/mpfr/2009-06/msg00000.html https://sympa.inria.fr/sympa/arc/mpfr/2009-06/msg00008.html https://sympa.inria.fr/sympa/arc/mpfr/2009-06/msg00010.html */ int mpfr_subnormalize (mpfr_ptr y, int old_inexact, mpfr_rnd_t rnd) { int sign; /* The subnormal exponent range is [ emin, emin + MPFR_PREC(y) - 2 ] */ if (MPFR_LIKELY (MPFR_IS_SINGULAR (y) || (MPFR_GET_EXP (y) >= __gmpfr_emin + (mpfr_exp_t) MPFR_PREC (y) - 1))) MPFR_RET (old_inexact); MPFR_SET_UNDERFLOW (); sign = MPFR_SIGN (y); /* We have to emulate one bit rounding if EXP(y) = emin */ if (MPFR_GET_EXP (y) == __gmpfr_emin) { /* If this is a power of 2, we don't need rounding. It handles cases when |y| = 0.1 * 2^emin */ if (mpfr_powerof2_raw (y)) MPFR_RET (old_inexact); /* We keep the same sign for y. Assuming Y is the real value and y the approximation and since y is not a power of 2: 0.5*2^emin < Y < 1*2^emin We also know the direction of the error thanks to ternary value. */ if (rnd == MPFR_RNDN || rnd == MPFR_RNDNA) { mp_limb_t *mant, rb ,sb; mp_size_t s; /* We need the rounding bit and the sticky bit. Read them and use the previous table to conclude. */ s = MPFR_LIMB_SIZE (y) - 1; mant = MPFR_MANT (y) + s; rb = *mant & (MPFR_LIMB_HIGHBIT >> 1); if (rb == 0) goto set_min; sb = *mant & ((MPFR_LIMB_HIGHBIT >> 1) - 1); while (sb == 0 && s-- != 0) sb = *--mant; if (sb != 0) goto set_min_p1; /* Rounding bit is 1 and sticky bit is 0. We need to examine old inexact flag to conclude. */ if ((old_inexact > 0 && sign > 0) || (old_inexact < 0 && sign < 0)) goto set_min; /* If inexact != 0, return 0.1*2^(emin+1). Otherwise, rounding bit = 1, sticky bit = 0 and inexact = 0 So we have 0.1100000000000000000000000*2^emin exactly. We return 0.1*2^(emin+1) according to the even-rounding rule on subnormals. Note the same holds for RNDNA. */ goto set_min_p1; } else if (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG (y))) { set_min: mpfr_setmin (y, __gmpfr_emin); MPFR_RET (-sign); } else { set_min_p1: /* Note: mpfr_setmin will abort if __gmpfr_emax == __gmpfr_emin. */ mpfr_setmin (y, __gmpfr_emin + 1); MPFR_RET (sign); } } else /* Hard case: It is more or less the same problem as mpfr_cache */ { mpfr_t dest; mpfr_prec_t q; mpfr_rnd_t rnd2; int inexact, inex2; MPFR_ASSERTD (MPFR_GET_EXP (y) > __gmpfr_emin); /* Compute the intermediary precision */ q = (mpfr_uexp_t) MPFR_GET_EXP (y) - __gmpfr_emin + 1; MPFR_ASSERTD (q >= MPFR_PREC_MIN && q < MPFR_PREC (y)); /* TODO: perform the rounding in place. */ mpfr_init2 (dest, q); /* Round y in dest */ MPFR_SET_EXP (dest, MPFR_GET_EXP (y)); MPFR_SET_SIGN (dest, sign); rnd2 = rnd == MPFR_RNDNA ? MPFR_RNDN : rnd; MPFR_RNDRAW_EVEN (inexact, dest, MPFR_MANT (y), MPFR_PREC (y), rnd2, sign, MPFR_SET_EXP (dest, MPFR_GET_EXP (dest) + 1)); if (MPFR_LIKELY (old_inexact != 0)) { if (MPFR_UNLIKELY (rnd2 == MPFR_RNDN && (inexact == MPFR_EVEN_INEX || inexact == -MPFR_EVEN_INEX))) { /* If both roundings are in the same direction, we have to go back in the other direction. For MPFR_RNDNA it is the same, since we are not exactly in the middle case (old_inexact != 0). */ if (SAME_SIGN (inexact, old_inexact)) { if (SAME_SIGN (inexact, MPFR_INT_SIGN (y))) mpfr_nexttozero (dest); else /* subnormal range, thus no overflow */ { mpfr_nexttoinf (dest); MPFR_ASSERTD(!MPFR_IS_INF (dest)); } inexact = -inexact; } } else if (MPFR_UNLIKELY (inexact == 0)) inexact = old_inexact; } else if (rnd == MPFR_RNDNA && (inexact == MPFR_EVEN_INEX || inexact == -MPFR_EVEN_INEX)) { /* We are in the middle case: since we used RNDN to round, we should round in the opposite direction when inexact has the opposite sign of y. */ if (!SAME_SIGN (inexact, MPFR_INT_SIGN (y))) { mpfr_nexttoinf (dest); MPFR_ASSERTD(!MPFR_IS_INF (dest)); inexact = -inexact; } } inex2 = mpfr_set (y, dest, rnd); MPFR_ASSERTN (inex2 == 0); MPFR_ASSERTN (MPFR_IS_PURE_FP (y)); mpfr_clear (dest); MPFR_RET (inexact); } }