"Fossies" - the Fresh Open Source Software Archive

Member "mpfr-4.0.2/src/subnormal.c" (7 Jan 2019, 6281 Bytes) of package /linux/misc/mpfr-4.0.2.tar.gz:


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 }