src/subnormal.c - master - githubmirror/mpfr

src/subnormal.c
210 строк8 KB
Vincent Lefevre
mpfr_subnormalize and underflow flag: final clean-up.
28 мар 2025, 14:06
28 мар 2025, 14:06fb00630
/* mpfr_subnormalize -- Subnormalize a floating point number
   emulating sub-normal numbers.

Copyright 2005-2025 Free Software Foundation, Inc.
Contributed by the Pascaline 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/>. */

#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 modes, 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 the first rounding (input of mpfr_subnormalize,
   with extended exponent range and default precision) is non-zero and is
   in the subnormal range.
   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

   Assuming that the user does not touch the underflow flag between
   the first rounding and the call to mpfr_subnormalize, this
   corresponds to "underflow after rounding", as explained below.
   Note that the underflow flag is set in all underflow cases,
   while in IEEE 754, under the default exception handling, the
   underflow flag is not set for the exact underflow cases.

   Explanation:
     * If the first rounding yields an underflow, then this will be a
       real underflow as the exponent range is larger than the one of
       the emulated system. The underflow flag is set during this first
       rounding.
     * If the first rounding does not yield an underflow, then this is
       like the computation were done with an unbounded exponent range,
       and mpfr_subnormalize sets the underflow flag only if its input
       is non-zero and is in the subnormal range, which corresponds to
       the definition of "underflow after rounding".

   Note that as documented, if the first rounding is an inexact zero,
   mpfr_subnormalize will not set the underflow flag. So, if the user
   has cleared the underflow flag after the first rounding, this will
   not correspond to "underflow after rounding".
*/

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 ].
     Note: if MPFR_PREC(y) is 1, the following condition is always true,
     which corresponds to the fact that there are no subnormals. */
  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.5 * 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.5*2^(emin+1).
             Otherwise, rounding bit = 1, sticky bit = 0 and inexact = 0
             So we have 0.1100000000000000000000000*2^emin exactly.
             We return 0.5*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);
    }
}