1 /* mpfr_get_float128 -- convert a multiple precision floating-point
2 number to a _Float128 number
3
4 Copyright 2012-2023 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramba 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 https://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 "mpfr-impl.h"
25
26 #ifdef MPFR_WANT_FLOAT128
27
28 /* generic code */
29 _Float128
mpfr_get_float128(mpfr_srcptr x,mpfr_rnd_t rnd_mode)30 mpfr_get_float128 (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
31 {
32
33 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
34 return (_Float128) mpfr_get_d (x, rnd_mode);
35 else /* now x is a normal non-zero number */
36 {
37 _Float128 r; /* result */
38 _Float128 m;
39 mpfr_exp_t e; /* exponent of x (before rounding) */
40 mpfr_exp_t sh; /* exponent shift, so that x/2^sh is in the double range */
41 const int emin = -16381;
42 const int esub = emin - IEEE_FLOAT128_MANT_DIG;
43 int sign;
44
45 sign = MPFR_SIGN (x);
46 e = MPFR_GET_EXP (x);
47
48 if (MPFR_UNLIKELY (e <= esub))
49 {
50 if (MPFR_IS_LIKE_RNDZ (rnd_mode, sign < 0) ||
51 (rnd_mode == MPFR_RNDN && (e < esub || mpfr_powerof2_raw (x))))
52 return sign < 0 ? -0.0 : 0.0;
53 r = 1.0;
54 sh = esub;
55 }
56 else
57 {
58 mpfr_t y;
59 mp_limb_t *yp;
60 int prec, i; /* small enough to fit in an int */
61 MPFR_SAVE_EXPO_DECL (expo);
62
63 MPFR_SAVE_EXPO_MARK (expo);
64
65 /* First round x to the target _Float128 precision, taking the
66 reduced precision of the subnormals into account, so that all
67 subsequent operations are exact (this avoids double rounding
68 problems). */
69 prec = e < emin ? e - esub : IEEE_FLOAT128_MANT_DIG;
70 MPFR_ASSERTD (prec >= MPFR_PREC_MIN);
71 mpfr_init2 (y, prec);
72
73 mpfr_set (y, x, rnd_mode);
74 sh = MPFR_GET_EXP (y);
75 MPFR_SET_EXP (y, 0);
76 MPFR_SET_POS (y);
77 yp = MPFR_MANT (y);
78
79 r = 0.0;
80 for (i = 0; i < MPFR_LIMB_SIZE (y); i++)
81 {
82 /* Note: MPFR_LIMB_MAX is avoided below as it might not
83 always work if GMP_NUMB_BITS > IEEE_FLOAT128_MANT_DIG.
84 MPFR_LIMB_HIGHBIT has the advantage to fit on 1 bit. */
85 r += yp[i];
86 r *= 1 / (2 * (_Float128) MPFR_LIMB_HIGHBIT);
87 }
88
89 mpfr_clear (y);
90
91 MPFR_SAVE_EXPO_FREE (expo);
92 }
93
94 /* we now have to multiply r by 2^sh */
95 MPFR_ASSERTD (r > 0);
96 if (sh != 0)
97 {
98 /* An overflow may occur (example: 0.5*2^1024) */
99 while (r < 1.0)
100 {
101 r += r;
102 sh--;
103 }
104
105 if (sh > 0)
106 m = 2.0;
107 else
108 {
109 m = 0.5;
110 sh = -sh;
111 }
112
113 for (;;)
114 {
115 if (sh % 2)
116 r = r * m;
117 sh >>= 1;
118 if (sh == 0)
119 break;
120 m = m * m;
121 }
122 }
123 if (sign < 0)
124 r = -r;
125 return r;
126 }
127 }
128
129 #endif /* MPFR_WANT_FLOAT128 */
130