Imported Upstream version 3.1.1
[platform/upstream/mpfr.git] / 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 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 }