1 /* Copyright (C) 1995-2023 Free Software Foundation, Inc.
2 This file is part of the GNU C Library.
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.
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.
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/>. */
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. */
31 __mpn_construct_long_double (mp_srcptr frac_ptr
, int expt
, int sign
)
33 union ibm_extended_long_double u
;
34 unsigned long lzcount
;
35 unsigned long long hi
, lo
;
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
;
44 #if BITS_PER_MP_LIMB == 32
45 /* The low order 53 bits (52 + hidden) go into the lower double */
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);
59 #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
62 if ((hi
& (1LL << 52)) == 0 && (hi
| lo
) != 0)
65 unsigned long long val
= hi
? hi
: lo
;
67 if (sizeof (val
) == sizeof (long))
68 lzcount
= __builtin_clzl (val
);
69 else if ((val
>> 32) != 0)
70 lzcount
= __builtin_clzl ((long) (val
>> 32));
72 lzcount
= __builtin_clzl ((long) val
) + 32;
74 lzcount
= lzcount
- (64 - 53);
76 lzcount
= lzcount
+ 53 - (64 - 53);
78 if (lzcount
> u
.d
[0].ieee
.exponent
)
80 lzcount
= u
.d
[0].ieee
.exponent
;
81 u
.d
[0].ieee
.exponent
= 0;
86 u
.d
[0].ieee
.exponent
-= (lzcount
- 1);
87 exponent2
-= (lzcount
- 1);
92 hi
= (hi
<< lzcount
) | (lo
>> (53 - lzcount
));
93 lo
= (lo
<< lzcount
) & ((1LL << 53) - 1);
97 hi
= lo
<< (lzcount
- 53);
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))
112 if ((hi
& (1LL << 53)) != 0)
115 u
.d
[0].ieee
.exponent
++;
116 if (u
.d
[0].ieee
.exponent
== IEEE754_DOUBLE_BIAS
+ DBL_MAX_EXP
)
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
;
125 u
.d
[1].ieee
.negative
= !sign
;
126 lo
= (1LL << 53) - lo
;
129 /* Normalize the low double. Shift the mantissa left until
130 the hidden bit is '1' and adjust the exponent accordingly. */
132 if (sizeof (lo
) == sizeof (long))
133 lzcount
= __builtin_clzl (lo
);
134 else if ((lo
>> 32) != 0)
135 lzcount
= __builtin_clzl ((long) (lo
>> 32));
137 lzcount
= __builtin_clzl ((long) lo
) + 32;
138 lzcount
= lzcount
- (64 - 53);
140 exponent2
-= lzcount
;
143 u
.d
[1].ieee
.exponent
= exponent2
;
144 else if (exponent2
> -53)
145 lo
>>= 1 - exponent2
;
150 u
.d
[1].ieee
.negative
= 0;
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;