1 /* Copyright (C) 1995-2023 Free Software Foundation, Inc.
2 This file is part of the GNU C Library.
3
4 The GNU C Library is free software; you can redistribute it and/or
5 modify it under the terms of the GNU Lesser General Public
6 License as published by the Free Software Foundation; either
7 version 2.1 of the License, or (at your option) any later version.
8
9 The GNU C Library is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 Lesser General Public License for more details.
13
14 You should have received a copy of the GNU Lesser General Public
15 License along with the GNU C Library; if not, see
16 <https://www.gnu.org/licenses/>. */
17
18 #include <ieee754.h>
19 #include <errno.h>
20 #include <float.h>
21 #include <math.h>
22
23 #include "gmp.h"
24 #include "gmp-impl.h"
25
26 /* Convert a multi-precision integer of the needed number of bits (106
27 for long double) and an integral power of two to a `long double' in
28 IBM extended format. */
29
30 long double
31 __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
32 {
33 union ibm_extended_long_double u;
34 unsigned long lzcount;
35 unsigned long long hi, lo;
36 int exponent2;
37
38 u.d[0].ieee.negative = sign;
39 u.d[1].ieee.negative = sign;
40 u.d[0].ieee.exponent = expt + IEEE754_DOUBLE_BIAS;
41 u.d[1].ieee.exponent = 0;
42 exponent2 = expt - 53 + IEEE754_DOUBLE_BIAS;
43
44 #if BITS_PER_MP_LIMB == 32
45 /* The low order 53 bits (52 + hidden) go into the lower double */
46 lo = frac_ptr[0];
47 lo |= (frac_ptr[1] & ((1LL << (53 - 32)) - 1)) << 32;
48 /* The high order 53 bits (52 + hidden) go into the upper double */
49 hi = (frac_ptr[1] >> (53 - 32)) & ((1 << 11) - 1);
50 hi |= ((unsigned long long) frac_ptr[2]) << 11;
51 hi |= ((unsigned long long) frac_ptr[3]) << (32 + 11);
52 #elif BITS_PER_MP_LIMB == 64
53 /* The low order 53 bits (52 + hidden) go into the lower double */
54 lo = frac_ptr[0] & (((mp_limb_t) 1 << 53) - 1);
55 /* The high order 53 bits (52 + hidden) go into the upper double */
56 hi = (frac_ptr[0] >> 53) & (((mp_limb_t) 1 << 11) - 1);
57 hi |= (frac_ptr[1] << 11);
58 #else
59 #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
60 #endif
61
62 if ((hi & (1LL << 52)) == 0 && (hi | lo) != 0)
63 {
64 /* denormal number */
65 unsigned long long val = hi ? hi : lo;
66
67 if (sizeof (val) == sizeof (long))
68 lzcount = __builtin_clzl (val);
69 else if ((val >> 32) != 0)
70 lzcount = __builtin_clzl ((long) (val >> 32));
71 else
72 lzcount = __builtin_clzl ((long) val) + 32;
73 if (hi)
74 lzcount = lzcount - (64 - 53);
75 else
76 lzcount = lzcount + 53 - (64 - 53);
77
78 if (lzcount > u.d[0].ieee.exponent)
79 {
80 lzcount = u.d[0].ieee.exponent;
81 u.d[0].ieee.exponent = 0;
82 exponent2 -= lzcount;
83 }
84 else
85 {
86 u.d[0].ieee.exponent -= (lzcount - 1);
87 exponent2 -= (lzcount - 1);
88 }
89
90 if (lzcount <= 53)
91 {
92 hi = (hi << lzcount) | (lo >> (53 - lzcount));
93 lo = (lo << lzcount) & ((1LL << 53) - 1);
94 }
95 else
96 {
97 hi = lo << (lzcount - 53);
98 lo = 0;
99 }
100 }
101
102 if (lo != 0)
103 {
104 /* hidden bit of low double controls rounding of the high double.
105 If hidden is '1' and either the explicit mantissa is non-zero
106 or hi is odd, then round up hi and adjust lo (2nd mantissa)
107 plus change the sign of the low double to compensate. */
108 if ((lo & (1LL << 52)) != 0
109 && ((hi & 1) != 0 || (lo & ((1LL << 52) - 1)) != 0))
110 {
111 hi++;
112 if ((hi & (1LL << 53)) != 0)
113 {
114 hi >>= 1;
115 u.d[0].ieee.exponent++;
116 if (u.d[0].ieee.exponent == IEEE754_DOUBLE_BIAS + DBL_MAX_EXP)
117 {
118 /* Overflow. The appropriate overflowed result must
119 be produced (if an infinity, that means the low
120 part must be zero). */
121 __set_errno (ERANGE);
122 return (sign ? -LDBL_MAX : LDBL_MAX) * LDBL_MAX;
123 }
124 }
125 u.d[1].ieee.negative = !sign;
126 lo = (1LL << 53) - lo;
127 }
128
129 /* Normalize the low double. Shift the mantissa left until
130 the hidden bit is '1' and adjust the exponent accordingly. */
131
132 if (sizeof (lo) == sizeof (long))
133 lzcount = __builtin_clzl (lo);
134 else if ((lo >> 32) != 0)
135 lzcount = __builtin_clzl ((long) (lo >> 32));
136 else
137 lzcount = __builtin_clzl ((long) lo) + 32;
138 lzcount = lzcount - (64 - 53);
139 lo <<= lzcount;
140 exponent2 -= lzcount;
141
142 if (exponent2 > 0)
143 u.d[1].ieee.exponent = exponent2;
144 else if (exponent2 > -53)
145 lo >>= 1 - exponent2;
146 else
147 lo = 0;
148 }
149 else
150 u.d[1].ieee.negative = 0;
151
152 u.d[1].ieee.mantissa1 = lo;
153 u.d[1].ieee.mantissa0 = lo >> 32;
154 u.d[0].ieee.mantissa1 = hi;
155 u.d[0].ieee.mantissa0 = hi >> 32;
156
157 return u.ld;
158 }