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   int inex;
 111   MPFR_SAVE_EXPO_DECL (expo);
 112
 113   MPFR_SAVE_EXPO_MARK (expo);
 114
 115   mpfr_init2 (tmp, MPFR_LDBL_MANT_DIG);
 116   inex = mpfr_set (tmp, x, rnd_mode);
 117
 118   mpfr_set_emin (-16382-63);
 119   mpfr_set_emax (16384);
 120   mpfr_subnormalize (tmp, mpfr_check_range (tmp, inex, rnd_mode), rnd_mode);
 121   mpfr_prec_round (tmp, 64, GMP_RNDZ); /* exact */
 122   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tmp)))
 123     ld.ld = (long double) mpfr_get_d (tmp, rnd_mode);
 124   else
 125     {
 126       mp_limb_t *tmpmant;
 127       mp_exp_t e, denorm;
 128
 129       tmpmant = MPFR_MANT (tmp);
 130       e = MPFR_GET_EXP (tmp);
 131       /* the smallest normal number is 2^(-16382), which is 0.5*2^(-16381)
 132          in MPFR, thus any exponent <= -16382 corresponds to a subnormal
 133          number */
 134       denorm = MPFR_UNLIKELY (e <= -16382) ? - e - 16382 + 1 : 0;
 135 #if BITS_PER_MP_LIMB >= 64
 136       ld.s.manl = (tmpmant[0] >> denorm);
 137       ld.s.manh = (tmpmant[0] >> denorm) >> 32;
 138 #elif BITS_PER_MP_LIMB == 32
 139       if (MPFR_LIKELY (denorm == 0))
 140         {
 141           ld.s.manl = tmpmant[0];
 142           ld.s.manh = tmpmant[1];
 143         }
 144       else if (denorm < 32)
 145         {
 146           ld.s.manl = (tmpmant[0] >> denorm) | (tmpmant[1] << (32 - denorm));
 147           ld.s.manh = tmpmant[1] >> denorm;
 148         }
 149       else /* 32 <= denorm <= 64 */
 150         {
 151           ld.s.manl = tmpmant[1] >> (denorm - 32);
 152           ld.s.manh = 0;
 153         }
 154 #else
 155 # error "BITS_PER_MP_LIMB must be 32 or >= 64"
 156       /* Other values have never been supported anyway. */
 157 #endif
 158       if (MPFR_LIKELY (denorm == 0))
 159         {
 160           ld.s.exph = (e + 0x3FFE) >> 8;
 161           ld.s.expl = (e + 0x3FFE);
 162         }
 163       else
 164         ld.s.exph = ld.s.expl = 0;
 165       ld.s.sign = MPFR_IS_NEG (x);
 166     }
 167
 168   mpfr_clear (tmp);
 169   MPFR_SAVE_EXPO_FREE (expo);
 170   return ld.ld;
 171 }
 172
 173 #endif
 174
 175 /* contributed by Damien Stehle */
 176 long double
 177 mpfr_get_ld_2exp (long *expptr, mpfr_srcptr src, mp_rnd_t rnd_mode)
 178 {
 179   long double ret;
 180   mp_exp_t exp;
 181   mpfr_t tmp;
 182
 183   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
 184     return (long double) mpfr_get_d_2exp (expptr, src, rnd_mode);
 185
 186   tmp[0] = *src;        /* Hack copy mpfr_t */
 187   MPFR_SET_EXP (tmp, 0);
 188   ret = mpfr_get_ld (tmp, rnd_mode);
 189
 190   if (MPFR_IS_PURE_FP(src))
 191     {
 192       exp = MPFR_GET_EXP (src);
 193
 194       /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */
 195       if (ret == 1.0)
 196         {
 197           ret = 0.5;
 198           exp ++;
 199         }
 200       else if (ret ==  -1.0)
 201         {
 202           ret = -0.5;
 203           exp ++;
 204         }
 205
 206       MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0)
 207                     || (ret <= -0.5 && ret > -1.0));
 208       MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX);
 209     }
 210   else
 211     exp = 0;
 212
 213   *expptr = exp;
 214   return ret;
 215 }