contrib/mpfr/src/get_ld.c

   1 /* mpfr_get_ld, mpfr_get_ld_2exp -- convert a multiple precision floating-point
   2                                     number to a machine long double
   3
   4 Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
   5 Contributed by the AriC and Caramel 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 http://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 <float.h>
  25
  26 #include "mpfr-impl.h"
  27
  28 #ifndef HAVE_LDOUBLE_IEEE_EXT_LITTLE
  29
  30 long double
  31 mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
  32 {
  33
  34   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
  35     return (long double) mpfr_get_d (x, rnd_mode);
  36   else /* now x is a normal non-zero number */
  37     {
  38       long double r; /* result */
  39       long double m;
  40       double s; /* part of result */
  41       mpfr_exp_t sh; /* exponent shift, so that x/2^sh is in the double range */
  42       mpfr_t y, z;
  43       int sign;
  44
  45       /* first round x to the target long double precision, so that
  46          all subsequent operations are exact (this avoids double rounding
  47          problems) */
  48       mpfr_init2 (y, MPFR_LDBL_MANT_DIG);
  49       mpfr_init2 (z, MPFR_LDBL_MANT_DIG);
  50       /* Note about the precision of z: even though IEEE_DBL_MANT_DIG is
  51          sufficient, z has been set to the same precision as y so that
  52          the mpfr_sub below calls mpfr_sub1sp, which is faster than the
  53          generic subtraction, even in this particular case (from tests
  54          done by Patrick Pelissier on a 64-bit Core2 Duo against r7285).
  55          But here there is an important cancellation in the subtraction.
  56          TODO: get more information about what has been tested. */
  57
  58       mpfr_set (y, x, rnd_mode);
  59       sh = MPFR_GET_EXP (y);
  60       sign = MPFR_SIGN (y);
  61       MPFR_SET_EXP (y, 0);
  62       MPFR_SET_POS (y);
  63
  64       r = 0.0;
  65       do {
  66         s = mpfr_get_d (y, MPFR_RNDN); /* high part of y */
  67         r += (long double) s;
  68         mpfr_set_d (z, s, MPFR_RNDN);  /* exact */
  69         mpfr_sub (y, y, z, MPFR_RNDN); /* exact */
  70       } while (!MPFR_IS_ZERO (y));
  71
  72       mpfr_clear (z);
  73       mpfr_clear (y);
  74
  75       /* we now have to multiply back by 2^sh */
  76       MPFR_ASSERTD (r > 0);
  77       if (sh != 0)
  78         {
  79           /* An overflow may occurs (example: 0.5*2^1024) */
  80           while (r < 1.0)
  81             {
  82               r += r;
  83               sh--;
  84             }
  85
  86           if (sh > 0)
  87             m = 2.0;
  88           else
  89             {
  90               m = 0.5;
  91               sh = -sh;
  92             }
  93
  94           for (;;)
  95             {
  96               if (sh % 2)
  97                 r = r * m;
  98               sh >>= 1;
  99               if (sh == 0)
 100                 break;
 101               m = m * m;
 102             }
 103         }
 104       if (sign < 0)
 105         r = -r;
 106       return r;
 107     }
 108 }
 109
 110 #else
 111
 112 long double
 113 mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
 114 {
 115   mpfr_long_double_t ld;
 116   mpfr_t tmp;
 117   int inex;
 118   MPFR_SAVE_EXPO_DECL (expo);
 119
 120   MPFR_SAVE_EXPO_MARK (expo);
 121
 122   mpfr_init2 (tmp, MPFR_LDBL_MANT_DIG);
 123   inex = mpfr_set (tmp, x, rnd_mode);
 124
 125   mpfr_set_emin (-16382-63);
 126   mpfr_set_emax (16384);
 127   mpfr_subnormalize (tmp, mpfr_check_range (tmp, inex, rnd_mode), rnd_mode);
 128   mpfr_prec_round (tmp, 64, MPFR_RNDZ); /* exact */
 129   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tmp)))
 130     ld.ld = (long double) mpfr_get_d (tmp, rnd_mode);
 131   else
 132     {
 133       mp_limb_t *tmpmant;
 134       mpfr_exp_t e, denorm;
 135
 136       tmpmant = MPFR_MANT (tmp);
 137       e = MPFR_GET_EXP (tmp);
 138       /* the smallest normal number is 2^(-16382), which is 0.5*2^(-16381)
 139          in MPFR, thus any exponent <= -16382 corresponds to a subnormal
 140          number */
 141       denorm = MPFR_UNLIKELY (e <= -16382) ? - e - 16382 + 1 : 0;
 142 #if GMP_NUMB_BITS >= 64
 143       ld.s.manl = (tmpmant[0] >> denorm);
 144       ld.s.manh = (tmpmant[0] >> denorm) >> 32;
 145 #elif GMP_NUMB_BITS == 32
 146       if (MPFR_LIKELY (denorm == 0))
 147         {
 148           ld.s.manl = tmpmant[0];
 149           ld.s.manh = tmpmant[1];
 150         }
 151       else if (denorm < 32)
 152         {
 153           ld.s.manl = (tmpmant[0] >> denorm) | (tmpmant[1] << (32 - denorm));
 154           ld.s.manh = tmpmant[1] >> denorm;
 155         }
 156       else /* 32 <= denorm <= 64 */
 157         {
 158           ld.s.manl = tmpmant[1] >> (denorm - 32);
 159           ld.s.manh = 0;
 160         }
 161 #else
 162 # error "GMP_NUMB_BITS must be 32 or >= 64"
 163       /* Other values have never been supported anyway. */
 164 #endif
 165       if (MPFR_LIKELY (denorm == 0))
 166         {
 167           ld.s.exph = (e + 0x3FFE) >> 8;
 168           ld.s.expl = (e + 0x3FFE);
 169         }
 170       else
 171         ld.s.exph = ld.s.expl = 0;
 172       ld.s.sign = MPFR_IS_NEG (x);
 173     }
 174
 175   mpfr_clear (tmp);
 176   MPFR_SAVE_EXPO_FREE (expo);
 177   return ld.ld;
 178 }
 179
 180 #endif
 181
 182 /* contributed by Damien Stehle */
 183 long double
 184 mpfr_get_ld_2exp (long *expptr, mpfr_srcptr src, mpfr_rnd_t rnd_mode)
 185 {
 186   long double ret;
 187   mpfr_exp_t exp;
 188   mpfr_t tmp;
 189
 190   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
 191     return (long double) mpfr_get_d_2exp (expptr, src, rnd_mode);
 192
 193   tmp[0] = *src;        /* Hack copy mpfr_t */
 194   MPFR_SET_EXP (tmp, 0);
 195   ret = mpfr_get_ld (tmp, rnd_mode);
 196
 197   if (MPFR_IS_PURE_FP(src))
 198     {
 199       exp = MPFR_GET_EXP (src);
 200
 201       /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */
 202       if (ret == 1.0)
 203         {
 204           ret = 0.5;
 205           exp ++;
 206         }
 207       else if (ret ==  -1.0)
 208         {
 209           ret = -0.5;
 210           exp ++;
 211         }
 212
 213       MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0)
 214                     || (ret <= -0.5 && ret > -1.0));
 215       MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX);
 216     }
 217   else
 218     exp = 0;
 219
 220   *expptr = exp;
 221   return ret;
 222 }