Revert not defining NO_HIDDEN on hurd
[glibc.git] / sysdeps / ieee754 / ldbl-128ibm / mpn2ldbl.c
blobc7ebd7ebd82e88a493163836b3987f3e4e6b5227
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/>. */
18 #include <ieee754.h>
19 #include <errno.h>
20 #include <float.h>
21 #include <math.h>
23 /* Need to set this when including gmp headers after system headers. */
24 #define HAVE_ALLOCA 1
26 #include "gmp.h"
27 #include "gmp-impl.h"
29 /* Convert a multi-precision integer of the needed number of bits (106
30 for long double) and an integral power of two to a `long double' in
31 IBM extended format. */
33 long double
34 __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
36 union ibm_extended_long_double u;
37 unsigned long lzcount;
38 unsigned long long hi, lo;
39 int exponent2;
41 u.d[0].ieee.negative = sign;
42 u.d[1].ieee.negative = sign;
43 u.d[0].ieee.exponent = expt + IEEE754_DOUBLE_BIAS;
44 u.d[1].ieee.exponent = 0;
45 exponent2 = expt - 53 + IEEE754_DOUBLE_BIAS;
47 #if BITS_PER_MP_LIMB == 32
48 /* The low order 53 bits (52 + hidden) go into the lower double */
49 lo = frac_ptr[0];
50 lo |= (frac_ptr[1] & ((1LL << (53 - 32)) - 1)) << 32;
51 /* The high order 53 bits (52 + hidden) go into the upper double */
52 hi = (frac_ptr[1] >> (53 - 32)) & ((1 << 11) - 1);
53 hi |= ((unsigned long long) frac_ptr[2]) << 11;
54 hi |= ((unsigned long long) frac_ptr[3]) << (32 + 11);
55 #elif BITS_PER_MP_LIMB == 64
56 /* The low order 53 bits (52 + hidden) go into the lower double */
57 lo = frac_ptr[0] & (((mp_limb_t) 1 << 53) - 1);
58 /* The high order 53 bits (52 + hidden) go into the upper double */
59 hi = (frac_ptr[0] >> 53) & (((mp_limb_t) 1 << 11) - 1);
60 hi |= (frac_ptr[1] << 11);
61 #else
62 #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
63 #endif
65 if ((hi & (1LL << 52)) == 0 && (hi | lo) != 0)
67 /* denormal number */
68 unsigned long long val = hi ? hi : lo;
70 if (sizeof (val) == sizeof (long))
71 lzcount = __builtin_clzl (val);
72 else if ((val >> 32) != 0)
73 lzcount = __builtin_clzl ((long) (val >> 32));
74 else
75 lzcount = __builtin_clzl ((long) val) + 32;
76 if (hi)
77 lzcount = lzcount - (64 - 53);
78 else
79 lzcount = lzcount + 53 - (64 - 53);
81 if (lzcount > u.d[0].ieee.exponent)
83 lzcount = u.d[0].ieee.exponent;
84 u.d[0].ieee.exponent = 0;
85 exponent2 -= lzcount;
87 else
89 u.d[0].ieee.exponent -= (lzcount - 1);
90 exponent2 -= (lzcount - 1);
93 if (lzcount <= 53)
95 hi = (hi << lzcount) | (lo >> (53 - lzcount));
96 lo = (lo << lzcount) & ((1LL << 53) - 1);
98 else
100 hi = lo << (lzcount - 53);
101 lo = 0;
105 if (lo != 0)
107 /* hidden bit of low double controls rounding of the high double.
108 If hidden is '1' and either the explicit mantissa is non-zero
109 or hi is odd, then round up hi and adjust lo (2nd mantissa)
110 plus change the sign of the low double to compensate. */
111 if ((lo & (1LL << 52)) != 0
112 && ((hi & 1) != 0 || (lo & ((1LL << 52) - 1)) != 0))
114 hi++;
115 if ((hi & (1LL << 53)) != 0)
117 hi >>= 1;
118 u.d[0].ieee.exponent++;
119 if (u.d[0].ieee.exponent == IEEE754_DOUBLE_BIAS + DBL_MAX_EXP)
121 /* Overflow. The appropriate overflowed result must
122 be produced (if an infinity, that means the low
123 part must be zero). */
124 __set_errno (ERANGE);
125 return (sign ? -LDBL_MAX : LDBL_MAX) * LDBL_MAX;
128 u.d[1].ieee.negative = !sign;
129 lo = (1LL << 53) - lo;
132 /* Normalize the low double. Shift the mantissa left until
133 the hidden bit is '1' and adjust the exponent accordingly. */
135 if (sizeof (lo) == sizeof (long))
136 lzcount = __builtin_clzl (lo);
137 else if ((lo >> 32) != 0)
138 lzcount = __builtin_clzl ((long) (lo >> 32));
139 else
140 lzcount = __builtin_clzl ((long) lo) + 32;
141 lzcount = lzcount - (64 - 53);
142 lo <<= lzcount;
143 exponent2 -= lzcount;
145 if (exponent2 > 0)
146 u.d[1].ieee.exponent = exponent2;
147 else if (exponent2 > -53)
148 lo >>= 1 - exponent2;
149 else
150 lo = 0;
152 else
153 u.d[1].ieee.negative = 0;
155 u.d[1].ieee.mantissa1 = lo;
156 u.d[1].ieee.mantissa0 = lo >> 32;
157 u.d[0].ieee.mantissa1 = hi;
158 u.d[0].ieee.mantissa0 = hi >> 32;
160 return u.ld;