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