contrib/mpfr/subnormal.c

   1 /* mpfr_subnormalize -- Subnormalize a floating point number
   2    emulating sub-normal numbers.
   3
   4 Copyright 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   5 Contributed by the Arenaire and Cacao 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 2.1 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.LIB.  If not, write to
  21 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  22 MA 02110-1301, USA. */
  23
  24 #include "mpfr-impl.h"
  25
  26 /* For GMP_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 inexact flags.
  37 */
  38
  39 int
  40 mpfr_subnormalize (mpfr_ptr y, int old_inexact, mp_rnd_t rnd)
  41 {
  42   int inexact = 0;
  43
  44   /* The subnormal exponent range are from:
  45       mpfr_emin to mpfr_emin + MPFR_PREC(y) - 1  */
  46   if (MPFR_LIKELY (MPFR_IS_SINGULAR (y)
  47                    || (MPFR_GET_EXP (y) >=
  48                        __gmpfr_emin + (mp_exp_t) MPFR_PREC (y) - 1)))
  49     inexact = old_inexact;
  50
  51   /* We have to emulate one bit rounding if EXP(y) = emin */
  52   else if (MPFR_GET_EXP (y) == __gmpfr_emin)
  53     {
  54       /* If this is a power of 2, we don't need rounding.
  55          It handles cases when rouding away and y=0.1*2^emin */
  56       if (mpfr_powerof2_raw (y))
  57         inexact = old_inexact;
  58       /* We keep the same sign for y.
  59          Assuming Y is the real value and y the approximation
  60          and since y is not a power of 2:  0.5*2^emin < Y < 1*2^emin
  61          We also know the direction of the error thanks to inexact flag */
  62       else if (rnd == GMP_RNDN)
  63         {
  64           mp_limb_t *mant, rb ,sb;
  65           mp_size_t s;
  66           /* We need the rounding bit and the sticky bit. Read them
  67              and use the previous table to conclude. */
  68           s = MPFR_LIMB_SIZE (y) - 1;
  69           mant = MPFR_MANT (y) + s;
  70           rb = *mant & (MPFR_LIMB_HIGHBIT >> 1);
  71           if (rb == 0)
  72             goto set_min;
  73           sb = *mant & ((MPFR_LIMB_HIGHBIT >> 1) - 1);
  74           while (sb == 0 && s-- != 0)
  75             sb = *--mant;
  76           if (sb != 0)
  77             goto set_min_p1;
  78           /* Rounding bit is 1 and sticky bit is 0.
  79              We need to examine old inexact flag to conclude. */
  80           if ((old_inexact > 0 && MPFR_SIGN (y) > 0) ||
  81               (old_inexact < 0 && MPFR_SIGN (y) < 0))
  82             goto set_min;
  83           /* If inexact != 0, return 0.1*2^(emin+1).
  84              Otherwise, rounding bit = 1, sticky bit = 0 and inexact = 0
  85              So we have 0.1100000000000000000000000*2^emin exactly.
  86              We return 0.1*2^(emin+1) according to the even-rounding
  87              rule on subnormals. */
  88           goto set_min_p1;
  89         }
  90       else if (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG (y)))
  91         {
  92         set_min:
  93           mpfr_setmin (y, __gmpfr_emin);
  94           inexact = -MPFR_SIGN (y);
  95         }
  96       else
  97         {
  98         set_min_p1:
  99           /* Note: mpfr_setmin will abort if __gmpfr_emax == __gmpfr_emin. */
 100           mpfr_setmin (y, __gmpfr_emin + 1);
 101           inexact = MPFR_SIGN (y);
 102         }
 103     }
 104
 105   else /* Hard case: It is more or less the same problem than mpfr_cache */
 106     {
 107       mpfr_t dest;
 108       mp_prec_t q;
 109       int sign;
 110
 111       /* Compute the intermediary precision */
 112       q = (mpfr_uexp_t) MPFR_GET_EXP (y) - __gmpfr_emin + 1;
 113       mpfr_init2 (dest, q);
 114       /* Round y in dest */
 115       sign = MPFR_SIGN (y);
 116       MPFR_SET_EXP (dest, MPFR_GET_EXP (y));
 117       MPFR_SET_SIGN (dest, sign);
 118       MPFR_RNDRAW_EVEN (inexact, dest,
 119                         MPFR_MANT (y), MPFR_PREC (y), rnd, sign,
 120                         MPFR_SET_EXP (dest, MPFR_GET_EXP (dest)+1));
 121       if (MPFR_LIKELY (old_inexact != 0))
 122         {
 123           if (MPFR_UNLIKELY(rnd == GMP_RNDN && (inexact == MPFR_EVEN_INEX
 124                                                || inexact == -MPFR_EVEN_INEX)))
 125             {
 126               /* if both roundings are in the same direction, we have to go
 127                  back in the other direction */
 128               if (SAME_SIGN (inexact, old_inexact))
 129                 {
 130                   if (SAME_SIGN (inexact, MPFR_INT_SIGN (y)))
 131                     mpfr_nexttozero (dest);
 132                   else
 133                     mpfr_nexttoinf (dest);
 134                   inexact = -inexact;
 135                 }
 136             }
 137           else if (MPFR_UNLIKELY (inexact == 0))
 138             inexact = old_inexact;
 139         }
 140       old_inexact = mpfr_set (y, dest, rnd);
 141       MPFR_ASSERTN (old_inexact == 0);
 142       MPFR_ASSERTN (MPFR_IS_PURE_FP (y));
 143       mpfr_clear (dest);
 144     }
 145   return inexact;
 146 }