1 /* Copyright (C) 1995-2015 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 <http://www.gnu.org/licenses/>. */
24 /* Convert a multi-precision integer of the needed number of bits (106
25 for long double) and an integral power of two to a `long double' in
26 IBM extended format. */
29 __mpn_construct_long_double (mp_srcptr frac_ptr
, int expt
, int sign
)
31 union ibm_extended_long_double u
;
32 unsigned long lzcount
;
33 unsigned long long hi
, lo
;
36 u
.d
[0].ieee
.negative
= sign
;
37 u
.d
[1].ieee
.negative
= sign
;
38 u
.d
[0].ieee
.exponent
= expt
+ IEEE754_DOUBLE_BIAS
;
39 u
.d
[1].ieee
.exponent
= 0;
40 exponent2
= expt
- 53 + IEEE754_DOUBLE_BIAS
;
42 #if BITS_PER_MP_LIMB == 32
43 /* The low order 53 bits (52 + hidden) go into the lower double */
45 lo
|= (frac_ptr
[1] & ((1LL << (53 - 32)) - 1)) << 32;
46 /* The high order 53 bits (52 + hidden) go into the upper double */
47 hi
= (frac_ptr
[1] >> (53 - 32)) & ((1 << 11) - 1);
48 hi
|= ((unsigned long long) frac_ptr
[2]) << 11;
49 hi
|= ((unsigned long long) frac_ptr
[3]) << (32 + 11);
50 #elif BITS_PER_MP_LIMB == 64
51 /* The low order 53 bits (52 + hidden) go into the lower double */
52 lo
= frac_ptr
[0] & (((mp_limb_t
) 1 << 53) - 1);
53 /* The high order 53 bits (52 + hidden) go into the upper double */
54 hi
= (frac_ptr
[0] >> 53) & (((mp_limb_t
) 1 << 11) - 1);
55 hi
|= (frac_ptr
[1] << 11);
57 #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
60 if ((hi
& (1LL << 52)) == 0 && (hi
| lo
) != 0)
63 unsigned long long val
= hi
? hi
: lo
;
65 if (sizeof (val
) == sizeof (long))
66 lzcount
= __builtin_clzl (val
);
67 else if ((val
>> 32) != 0)
68 lzcount
= __builtin_clzl ((long) (val
>> 32));
70 lzcount
= __builtin_clzl ((long) val
) + 32;
72 lzcount
= lzcount
- (64 - 53);
74 lzcount
= lzcount
+ 53 - (64 - 53);
76 if (lzcount
> u
.d
[0].ieee
.exponent
)
78 lzcount
= u
.d
[0].ieee
.exponent
;
79 u
.d
[0].ieee
.exponent
= 0;
84 u
.d
[0].ieee
.exponent
-= (lzcount
- 1);
85 exponent2
-= (lzcount
- 1);
90 hi
= (hi
<< lzcount
) | (lo
>> (53 - lzcount
));
91 lo
= (lo
<< lzcount
) & ((1LL << 53) - 1);
95 hi
= lo
<< (lzcount
- 53);
102 /* hidden bit of low double controls rounding of the high double.
103 If hidden is '1' and either the explicit mantissa is non-zero
104 or hi is odd, then round up hi and adjust lo (2nd mantissa)
105 plus change the sign of the low double to compensate. */
106 if ((lo
& (1LL << 52)) != 0
107 && ((hi
& 1) != 0 || (lo
& ((1LL << 52) - 1)) != 0))
110 if ((hi
& (1LL << 53)) != 0)
113 u
.d
[0].ieee
.exponent
++;
115 u
.d
[1].ieee
.negative
= !sign
;
116 lo
= (1LL << 53) - lo
;
119 /* Normalize the low double. Shift the mantissa left until
120 the hidden bit is '1' and adjust the exponent accordingly. */
122 if (sizeof (lo
) == sizeof (long))
123 lzcount
= __builtin_clzl (lo
);
124 else if ((lo
>> 32) != 0)
125 lzcount
= __builtin_clzl ((long) (lo
>> 32));
127 lzcount
= __builtin_clzl ((long) lo
) + 32;
128 lzcount
= lzcount
- (64 - 53);
130 exponent2
-= lzcount
;
133 u
.d
[1].ieee
.exponent
= exponent2
;
134 else if (exponent2
> -53)
135 lo
>>= 1 - exponent2
;
140 u
.d
[1].ieee
.negative
= 0;
142 u
.d
[1].ieee
.mantissa1
= lo
;
143 u
.d
[1].ieee
.mantissa0
= lo
>> 32;
144 u
.d
[0].ieee
.mantissa1
= hi
;
145 u
.d
[0].ieee
.mantissa0
= hi
>> 32;