1 /* Manipulation of the bit representation of 'long double' quantities.
2 Copyright (C) 2006-2023 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
16 License along with the GNU C Library; if not, see
17 <https://www.gnu.org/licenses/>. */
18
19 #ifndef _MATH_LDBL_H_
20 #define _MATH_LDBL_H_ 1
21
22 #include <ieee754.h>
23 #include <stdint.h>
24
25 /* To suit our callers we return *hi64 and *lo64 as if they came from
26 an ieee854 112 bit mantissa, that is, 48 bits in *hi64 (plus one
27 implicit bit) and 64 bits in *lo64. */
28
29 static inline void
30 ldbl_extract_mantissa (int64_t *hi64, uint64_t *lo64, int *exp, long double x)
31 {
32 /* We have 105 bits of mantissa plus one implicit digit. Since
33 106 bits are representable we use the first implicit digit for
34 the number before the decimal point and the second implicit bit
35 as bit 53 of the mantissa. */
36 uint64_t hi, lo;
37 union ibm_extended_long_double u;
38
39 u.ld = x;
40 *exp = u.d[0].ieee.exponent - IEEE754_DOUBLE_BIAS;
41
42 lo = ((uint64_t) u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1;
43 hi = ((uint64_t) u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1;
44
45 if (u.d[0].ieee.exponent != 0)
46 {
47 int ediff;
48
49 /* If not a denormal or zero then we have an implicit 53rd bit. */
50 hi |= (uint64_t) 1 << 52;
51
52 if (u.d[1].ieee.exponent != 0)
53 lo |= (uint64_t) 1 << 52;
54 else
55 /* A denormal is to be interpreted as having a biased exponent
56 of 1. */
57 lo = lo << 1;
58
59 /* We are going to shift 4 bits out of hi later, because we only
60 want 48 bits in *hi64. That means we want 60 bits in lo, but
61 we currently only have 53. Shift the value up. */
62 lo = lo << 7;
63
64 /* The lower double is normalized separately from the upper.
65 We may need to adjust the lower mantissa to reflect this.
66 The difference between the exponents can be larger than 53
67 when the low double is much less than 1ULP of the upper
68 (in which case there are significant bits, all 0's or all
69 1's, between the two significands). The difference between
70 the exponents can be less than 53 when the upper double
71 exponent is nearing its minimum value (in which case the low
72 double is denormal ie. has an exponent of zero). */
73 ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53;
74 if (ediff > 0)
75 {
76 if (ediff < 64)
77 lo = lo >> ediff;
78 else
79 lo = 0;
80 }
81 else if (ediff < 0)
82 lo = lo << -ediff;
83
84 if (u.d[0].ieee.negative != u.d[1].ieee.negative
85 && lo != 0)
86 {
87 hi--;
88 lo = ((uint64_t) 1 << 60) - lo;
89 if (hi < (uint64_t) 1 << 52)
90 {
91 /* We have a borrow from the hidden bit, so shift left 1. */
92 hi = (hi << 1) | (lo >> 59);
93 lo = (((uint64_t) 1 << 60) - 1) & (lo << 1);
94 *exp = *exp - 1;
95 }
96 }
97 }
98 else
99 /* If the larger magnitude double is denormal then the smaller
100 one must be zero. */
101 hi = hi << 1;
102
103 *lo64 = (hi << 60) | lo;
104 *hi64 = hi >> 4;
105 }
106
107 static inline long double
108 ldbl_insert_mantissa (int sign, int exp, int64_t hi64, uint64_t lo64)
109 {
110 union ibm_extended_long_double u;
111 int expnt2;
112 uint64_t hi, lo;
113
114 u.d[0].ieee.negative = sign;
115 u.d[1].ieee.negative = sign;
116 u.d[0].ieee.exponent = exp + IEEE754_DOUBLE_BIAS;
117 u.d[1].ieee.exponent = 0;
118 expnt2 = exp - 53 + IEEE754_DOUBLE_BIAS;
119
120 /* Expect 113 bits (112 bits + hidden) right justified in two longs.
121 The low order 53 bits (52 + hidden) go into the lower double */
122 lo = (lo64 >> 7) & (((uint64_t) 1 << 53) - 1);
123 /* The high order 53 bits (52 + hidden) go into the upper double */
124 hi = lo64 >> 60;
125 hi |= hi64 << 4;
126
127 if (lo != 0)
128 {
129 int lzcount;
130
131 /* hidden bit of low double controls rounding of the high double.
132 If hidden is '1' and either the explicit mantissa is non-zero
133 or hi is odd, then round up hi and adjust lo (2nd mantissa)
134 plus change the sign of the low double to compensate. */
135 if ((lo & ((uint64_t) 1 << 52)) != 0
136 && ((hi & 1) != 0 || (lo & (((uint64_t) 1 << 52) - 1)) != 0))
137 {
138 hi++;
139 if ((hi & ((uint64_t) 1 << 53)) != 0)
140 {
141 hi = hi >> 1;
142 u.d[0].ieee.exponent++;
143 }
144 u.d[1].ieee.negative = !sign;
145 lo = ((uint64_t) 1 << 53) - lo;
146 }
147
148 /* Normalize the low double. Shift the mantissa left until
149 the hidden bit is '1' and adjust the exponent accordingly. */
150
151 if (sizeof (lo) == sizeof (long))
152 lzcount = __builtin_clzl (lo);
153 else if ((lo >> 32) != 0)
154 lzcount = __builtin_clzl ((long) (lo >> 32));
155 else
156 lzcount = __builtin_clzl ((long) lo) + 32;
157 lzcount = lzcount - (64 - 53);
158 lo <<= lzcount;
159 expnt2 -= lzcount;
160
161 if (expnt2 >= 1)
162 /* Not denormal. */
163 u.d[1].ieee.exponent = expnt2;
164 else
165 {
166 /* Is denormal. Note that biased exponent of 0 is treated
167 as if it was 1, hence the extra shift. */
168 if (expnt2 > -53)
169 lo >>= 1 - expnt2;
170 else
171 lo = 0;
172 }
173 }
174 else
175 u.d[1].ieee.negative = 0;
176
177 u.d[1].ieee.mantissa1 = lo;
178 u.d[1].ieee.mantissa0 = lo >> 32;
179 u.d[0].ieee.mantissa1 = hi;
180 u.d[0].ieee.mantissa0 = hi >> 32;
181 return u.ld;
182 }
183
184 /* Handy utility functions to pack/unpack/cononicalize and find the nearbyint
185 of long double implemented as double double. */
186 static inline long double
187 default_ldbl_pack (double a, double aa)
188 {
189 union ibm_extended_long_double u;
190 u.d[0].d = a;
191 u.d[1].d = aa;
192 return u.ld;
193 }
194
195 static inline void
196 default_ldbl_unpack (long double l, double *a, double *aa)
197 {
198 union ibm_extended_long_double u;
199 u.ld = l;
200 *a = u.d[0].d;
201 *aa = u.d[1].d;
202 }
203
204 #ifndef ldbl_pack
205 # define ldbl_pack default_ldbl_pack
206 #endif
207 #ifndef ldbl_unpack
208 # define ldbl_unpack default_ldbl_unpack
209 #endif
210
211 /* Extract high double. */
212 #define ldbl_high(x) ((double) x)
213
214 /* Convert a finite long double to canonical form.
215 Does not handle +/-Inf properly. */
216 static inline void
217 ldbl_canonicalize (double *a, double *aa)
218 {
219 double xh, xl;
220
221 xh = *a + *aa;
222 xl = (*a - xh) + *aa;
223 *a = xh;
224 *aa = xl;
225 }
226
227 /* Simple inline nearbyint (double) function.
228 Only works in the default rounding mode
229 but is useful in long double rounding functions. */
230 static inline double
231 ldbl_nearbyint (double a)
232 {
233 double two52 = 0x1p52;
234
235 if (__glibc_likely ((__builtin_fabs (a) < two52)))
236 {
237 if (__glibc_likely ((a > 0.0)))
238 {
239 a += two52;
240 a -= two52;
241 }
242 else if (__glibc_likely ((a < 0.0)))
243 {
244 a = two52 - a;
245 a = -(a - two52);
246 }
247 }
248 return a;
249 }
250
251 /* Canonicalize a result from an integer rounding function, in any
252 rounding mode. *A and *AA are finite and integers, with *A being
253 nonzero; if the result is not already canonical, *AA is plus or
254 minus a power of 2 that does not exceed the least set bit in
255 *A. */
256 static inline void
257 ldbl_canonicalize_int (double *a, double *aa)
258 {
259 /* Previously we used EXTRACT_WORDS64 from math_private.h, but in order
260 to avoid including internal headers we duplicate that code here. */
261 uint64_t ax, aax;
262 union { double value; uint64_t word; } extractor;
263 extractor.value = *a;
264 ax = extractor.word;
265 extractor.value = *aa;
266 aax = extractor.word;
267
268 int expdiff = ((ax >> 52) & 0x7ff) - ((aax >> 52) & 0x7ff);
269 if (expdiff <= 53)
270 {
271 if (expdiff == 53)
272 {
273 /* Half way between two double values; noncanonical iff the
274 low bit of A's mantissa is 1. */
275 if ((ax & 1) != 0)
276 {
277 *a += 2 * *aa;
278 *aa = -*aa;
279 }
280 }
281 else
282 {
283 /* The sum can be represented in a single double. */
284 *a += *aa;
285 *aa = 0;
286 }
287 }
288 }
289
290 #endif /* math_ldbl.h */