Warning: This file is not a C or C++ file. It does not have highlighting.
| 1 | /* Manipulation of the bit representation of 'long double' quantities. |
|---|---|
| 2 | Copyright (C) 2006-2022 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 */ |
| 291 |
Warning: This file is not a C or C++ file. It does not have highlighting.