1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
|
#ifndef _MATH_PRIVATE_H_
#error "Never use <math_ldbl.h> directly; include <math_private.h> instead."
#endif
#include <sysdeps/ieee754/ldbl-128/math_ldbl.h>
#include <ieee754.h>
#include <stdint.h>
/* To suit our callers we return *hi64 and *lo64 as if they came from
an ieee854 112 bit mantissa, that is, 48 bits in *hi64 (plus one
implicit bit) and 64 bits in *lo64. */
static inline void
ldbl_extract_mantissa (int64_t *hi64, uint64_t *lo64, int *exp, long double x)
{
/* We have 105 bits of mantissa plus one implicit digit. Since
106 bits are representable we use the first implicit digit for
the number before the decimal point and the second implicit bit
as bit 53 of the mantissa. */
uint64_t hi, lo;
union ibm_extended_long_double u;
u.ld = x;
*exp = u.d[0].ieee.exponent - IEEE754_DOUBLE_BIAS;
lo = ((uint64_t) u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1;
hi = ((uint64_t) u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1;
if (u.d[0].ieee.exponent != 0)
{
int ediff;
/* If not a denormal or zero then we have an implicit 53rd bit. */
hi |= (uint64_t) 1 << 52;
if (u.d[1].ieee.exponent != 0)
lo |= (uint64_t) 1 << 52;
else
/* A denormal is to be interpreted as having a biased exponent
of 1. */
lo = lo << 1;
/* We are going to shift 4 bits out of hi later, because we only
want 48 bits in *hi64. That means we want 60 bits in lo, but
we currently only have 53. Shift the value up. */
lo = lo << 7;
/* The lower double is normalized separately from the upper.
We may need to adjust the lower mantissa to reflect this.
The difference between the exponents can be larger than 53
when the low double is much less than 1ULP of the upper
(in which case there are significant bits, all 0's or all
1's, between the two significands). The difference between
the exponents can be less than 53 when the upper double
exponent is nearing its minimum value (in which case the low
double is denormal ie. has an exponent of zero). */
ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53;
if (ediff > 0)
{
if (ediff < 64)
lo = lo >> ediff;
else
lo = 0;
}
else if (ediff < 0)
lo = lo << -ediff;
if (u.d[0].ieee.negative != u.d[1].ieee.negative
&& lo != 0)
{
hi--;
lo = ((uint64_t) 1 << 60) - lo;
if (hi < (uint64_t) 1 << 52)
{
/* We have a borrow from the hidden bit, so shift left 1. */
hi = (hi << 1) | (lo >> 59);
lo = (((uint64_t) 1 << 60) - 1) & (lo << 1);
*exp = *exp - 1;
}
}
}
else
/* If the larger magnitude double is denormal then the smaller
one must be zero. */
hi = hi << 1;
*lo64 = (hi << 60) | lo;
*hi64 = hi >> 4;
}
static inline long double
ldbl_insert_mantissa (int sign, int exp, int64_t hi64, uint64_t lo64)
{
union ibm_extended_long_double u;
int expnt2;
uint64_t hi, lo;
u.d[0].ieee.negative = sign;
u.d[1].ieee.negative = sign;
u.d[0].ieee.exponent = exp + IEEE754_DOUBLE_BIAS;
u.d[1].ieee.exponent = 0;
expnt2 = exp - 53 + IEEE754_DOUBLE_BIAS;
/* Expect 113 bits (112 bits + hidden) right justified in two longs.
The low order 53 bits (52 + hidden) go into the lower double */
lo = (lo64 >> 7) & (((uint64_t) 1 << 53) - 1);
/* The high order 53 bits (52 + hidden) go into the upper double */
hi = lo64 >> 60;
hi |= hi64 << 4;
if (lo != 0)
{
int lzcount;
/* hidden bit of low double controls rounding of the high double.
If hidden is '1' and either the explicit mantissa is non-zero
or hi is odd, then round up hi and adjust lo (2nd mantissa)
plus change the sign of the low double to compensate. */
if ((lo & ((uint64_t) 1 << 52)) != 0
&& ((hi & 1) != 0 || (lo & (((uint64_t) 1 << 52) - 1)) != 0))
{
hi++;
if ((hi & ((uint64_t) 1 << 53)) != 0)
{
hi = hi >> 1;
u.d[0].ieee.exponent++;
}
u.d[1].ieee.negative = !sign;
lo = ((uint64_t) 1 << 53) - lo;
}
/* Normalize the low double. Shift the mantissa left until
the hidden bit is '1' and adjust the exponent accordingly. */
if (sizeof (lo) == sizeof (long))
lzcount = __builtin_clzl (lo);
else if ((lo >> 32) != 0)
lzcount = __builtin_clzl ((long) (lo >> 32));
else
lzcount = __builtin_clzl ((long) lo) + 32;
lzcount = lzcount - (64 - 53);
lo <<= lzcount;
expnt2 -= lzcount;
if (expnt2 >= 1)
/* Not denormal. */
u.d[1].ieee.exponent = expnt2;
else
{
/* Is denormal. Note that biased exponent of 0 is treated
as if it was 1, hence the extra shift. */
if (expnt2 > -53)
lo >>= 1 - expnt2;
else
lo = 0;
}
}
else
u.d[1].ieee.negative = 0;
u.d[1].ieee.mantissa1 = lo;
u.d[1].ieee.mantissa0 = lo >> 32;
u.d[0].ieee.mantissa1 = hi;
u.d[0].ieee.mantissa0 = hi >> 32;
return u.ld;
}
/* Handy utility functions to pack/unpack/cononicalize and find the nearbyint
of long double implemented as double double. */
static inline long double
default_ldbl_pack (double a, double aa)
{
union ibm_extended_long_double u;
u.d[0].d = a;
u.d[1].d = aa;
return u.ld;
}
static inline void
default_ldbl_unpack (long double l, double *a, double *aa)
{
union ibm_extended_long_double u;
u.ld = l;
*a = u.d[0].d;
*aa = u.d[1].d;
}
#ifndef ldbl_pack
# define ldbl_pack default_ldbl_pack
#endif
#ifndef ldbl_unpack
# define ldbl_unpack default_ldbl_unpack
#endif
/* Extract high double. */
#define ldbl_high(x) ((double) x)
/* Convert a finite long double to canonical form.
Does not handle +/-Inf properly. */
static inline void
ldbl_canonicalize (double *a, double *aa)
{
double xh, xl;
xh = *a + *aa;
xl = (*a - xh) + *aa;
*a = xh;
*aa = xl;
}
/* Simple inline nearbyint (double) function.
Only works in the default rounding mode
but is useful in long double rounding functions. */
static inline double
ldbl_nearbyint (double a)
{
double two52 = 0x1p52;
if (__builtin_expect ((__builtin_fabs (a) < two52), 1))
{
if (__builtin_expect ((a > 0.0), 1))
{
a += two52;
a -= two52;
}
else if (__builtin_expect ((a < 0.0), 1))
{
a = two52 - a;
a = -(a - two52);
}
}
return a;
}
|
