contrib/mpfr/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 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 <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, mp_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       mp_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, IEEE_DBL_MANT_DIG);
  50
  51       mpfr_set (y, x, rnd_mode);
  52       sh = MPFR_GET_EXP (y);
  53       sign = MPFR_SIGN (y);
  54       MPFR_SET_EXP (y, 0);
  55       MPFR_SET_POS (y);
  56
  57       r = 0.0;
  58       do {
  59         s = mpfr_get_d (y, GMP_RNDN); /* high part of y */
  60         r += (long double) s;
  61         mpfr_set_d (z, s, GMP_RNDN);  /* exact */
  62         mpfr_sub (y, y, z, GMP_RNDN); /* exact */
  63       } while (!MPFR_IS_ZERO (y));
  64
  65       mpfr_clear (z);
  66       mpfr_clear (y);
  67
  68       /* we now have to multiply back by 2^sh */
  69       MPFR_ASSERTD (r > 0);
  70       if (sh != 0)
  71         {
  72           /* An overflow may occurs (example: 0.5*2^1024) */
  73           while (r < 1.0)
  74             {
  75               r += r;
  76               sh--;
  77             }
  78
  79           if (sh > 0)
  80             m = 2.0;
  81           else
  82             {
  83               m = 0.5;
  84               sh = -sh;
  85             }
  86
  87           for (;;)
  88             {
  89               if (sh % 2)
  90                 r = r * m;
  91               sh >>= 1;
  92               if (sh == 0)
  93                 break;
  94               m = m * m;
  95             }
  96         }
  97       if (sign < 0)
  98         r = -r;
  99       return r;
 100     }
 101 }
 102
 103 #else
 104
 105 long double
 106 mpfr_get_ld (mpfr_srcptr x, mp_rnd_t rnd_mode)
 107 {
 108   mpfr_long_double_t ld;
 109   mpfr_t tmp;
 110   MPFR_SAVE_EXPO_DECL (expo);
 111
 112   MPFR_SAVE_EXPO_MARK (expo);
 113   mpfr_set_emin (-16382-63);
 114   mpfr_set_emax (16383);
 115
 116   mpfr_init2 (tmp, MPFR_LDBL_MANT_DIG);
 117   mpfr_subnormalize(tmp, mpfr_set (tmp, x, rnd_mode), rnd_mode);
 118   mpfr_prec_round (tmp, 64, GMP_RNDZ); /* exact */
 119   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tmp)))
 120     ld.ld = (long double) mpfr_get_d (tmp, rnd_mode);
 121   else
 122     {
 123       mp_limb_t *tmpmant;
 124       mp_exp_t e, denorm;
 125
 126       tmpmant = MPFR_MANT (tmp);
 127       e = MPFR_GET_EXP (tmp);
 128       denorm = MPFR_UNLIKELY (e < -16382) ? - e - 16382 + 1 : 0;
 129 #if BITS_PER_MP_LIMB >= 64
 130       ld.s.manl = (tmpmant[0] >> denorm);
 131       ld.s.manh = (tmpmant[0] >> denorm) >> 32;
 132 #elif BITS_PER_MP_LIMB == 32
 133       if (MPFR_LIKELY (denorm == 0))
 134         {
 135           ld.s.manl = tmpmant[0];
 136           ld.s.manh = tmpmant[1];
 137         }
 138       else if (denorm < 32)
 139         {
 140           ld.s.manl = (tmpmant[0] >> denorm) | (tmpmant[1] << (32 - denorm));
 141           ld.s.manh = tmpmant[1] >> denorm;
 142         }
 143       else /* 32 <= denorm <= 64 */
 144         {
 145           ld.s.manl = tmpmant[1] >> (denorm - 32);
 146           ld.s.manh = 0;
 147         }
 148 #else
 149 # error "BITS_PER_MP_LIMB must be 32 or >= 64"
 150       /* Other values have never been supported anyway. */
 151 #endif
 152       if (MPFR_LIKELY (denorm == 0))
 153         {
 154           ld.s.exph = (e + 0x3FFE) >> 8;
 155           ld.s.expl = (e + 0x3FFE);
 156         }
 157       else
 158         ld.s.exph = ld.s.expl = 0;
 159       ld.s.sign = MPFR_IS_NEG (x);
 160     }
 161
 162   mpfr_clear (tmp);
 163   MPFR_SAVE_EXPO_FREE (expo);
 164   return ld.ld;
 165 }
 166
 167 #endif
 168
 169 /* contributed by Damien Stehle */
 170 long double
 171 mpfr_get_ld_2exp (long *expptr, mpfr_srcptr src, mp_rnd_t rnd_mode)
 172 {
 173   long double ret;
 174   mp_exp_t exp;
 175   mpfr_t tmp;
 176
 177   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
 178     return (long double) mpfr_get_d_2exp (expptr, src, rnd_mode);
 179
 180   tmp[0] = *src;        /* Hack copy mpfr_t */
 181   MPFR_SET_EXP (tmp, 0);
 182   ret = mpfr_get_ld (tmp, rnd_mode);
 183
 184   if (MPFR_IS_PURE_FP(src))
 185     {
 186       exp = MPFR_GET_EXP (src);
 187
 188       /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */
 189       if (ret == 1.0)
 190         {
 191           ret = 0.5;
 192           exp ++;
 193         }
 194       else if (ret ==  -1.0)
 195         {
 196           ret = -0.5;
 197           exp ++;
 198         }
 199
 200       MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0)
 201                     || (ret <= -0.5 && ret > -1.0));
 202       MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX);
 203     }
 204   else
 205     exp = 0;
 206
 207   *expptr = exp;
 208   return ret;
 209 }