| 1 | /* |
| 2 | * QEMU float support |
| 3 | * |
| 4 | * The code in this source file is derived from release 2a of the SoftFloat |
| 5 | * IEC/IEEE Floating-point Arithmetic Package. Those parts of the code (and |
| 6 | * some later contributions) are provided under that license, as detailed below. |
| 7 | * It has subsequently been modified by contributors to the QEMU Project, |
| 8 | * so some portions are provided under: |
| 9 | * the SoftFloat-2a license |
| 10 | * the BSD license |
| 11 | * GPL-v2-or-later |
| 12 | * |
| 13 | * Any future contributions to this file after December 1st 2014 will be |
| 14 | * taken to be licensed under the Softfloat-2a license unless specifically |
| 15 | * indicated otherwise. |
| 16 | */ |
| 17 | |
| 18 | FloatPartsN partsN(return_nan)(const FloatPartsN *a, float_status *s) |
| 19 | { |
| 20 | switch (a->cls) { |
| 21 | case float_class_snan: |
| 22 | float_raise(float_flag_invalid | float_flag_invalid_snan, s); |
| 23 | if (get_default_nan_mode(s)) { |
| 24 | return partsN(default_nan)(s); |
| 25 | } else { |
| 26 | return partsN(silence_nan)(a, s); |
| 27 | } |
| 28 | break; |
| 29 | case float_class_qnan: |
| 30 | if (get_default_nan_mode(s)) { |
| 31 | return partsN(default_nan)(s); |
| 32 | } |
| 33 | break; |
| 34 | default: |
| 35 | g_assert_not_reached(); |
| 36 | } |
| 37 | return *a; |
| 38 | } |
| 39 | |
| 40 | FloatPartsN partsN(pick_nan)(const FloatPartsN *a, const FloatPartsN *b, |
| 41 | float_status *s) |
| 42 | { |
| 43 | bool have_snan = false; |
| 44 | const FloatPartsN *ret; |
| 45 | int cmp; |
| 46 | |
| 47 | if (is_snan(a->cls) || is_snan(b->cls)) { |
| 48 | float_raise(float_flag_invalid | float_flag_invalid_snan, s); |
| 49 | have_snan = true; |
| 50 | } |
| 51 | |
| 52 | if (get_default_nan_mode(s)) { |
| 53 | return partsN(default_nan)(s); |
| 54 | } |
| 55 | |
| 56 | switch (get_float_2nan_prop_rule(s)) { |
| 57 | case float_2nan_prop_s_ab: |
| 58 | if (have_snan) { |
| 59 | ret = is_snan(a->cls) ? a : b; |
| 60 | break; |
| 61 | } |
| 62 | /* fall through */ |
| 63 | case float_2nan_prop_ab: |
| 64 | ret = is_nan(a->cls) ? a : b; |
| 65 | break; |
| 66 | case float_2nan_prop_s_ba: |
| 67 | if (have_snan) { |
| 68 | ret = is_snan(b->cls) ? b : a; |
| 69 | break; |
| 70 | } |
| 71 | /* fall through */ |
| 72 | case float_2nan_prop_ba: |
| 73 | ret = is_nan(b->cls) ? b : a; |
| 74 | break; |
| 75 | case float_2nan_prop_x87: |
| 76 | /* |
| 77 | * This implements x87 NaN propagation rules: |
| 78 | * SNaN + QNaN => return the QNaN |
| 79 | * two SNaNs => return the one with the larger significand, silenced |
| 80 | * two QNaNs => return the one with the larger significand |
| 81 | * SNaN and a non-NaN => return the SNaN, silenced |
| 82 | * QNaN and a non-NaN => return the QNaN |
| 83 | * |
| 84 | * If we get down to comparing significands and they are the same, |
| 85 | * return the NaN with the positive sign bit (if any). |
| 86 | */ |
| 87 | if (is_snan(a->cls)) { |
| 88 | if (!is_snan(b->cls)) { |
| 89 | ret = is_qnan(b->cls) ? b : a; |
| 90 | break; |
| 91 | } |
| 92 | } else if (is_qnan(a->cls)) { |
| 93 | if (is_snan(b->cls) || !is_qnan(b->cls)) { |
| 94 | ret = a; |
| 95 | break; |
| 96 | } |
| 97 | } else { |
| 98 | ret = b; |
| 99 | break; |
| 100 | } |
| 101 | cmp = fracN(cmp)(a, b); |
| 102 | if (cmp == 0) { |
| 103 | cmp = a->sign < b->sign; |
| 104 | } |
| 105 | ret = cmp > 0 ? a : b; |
| 106 | break; |
| 107 | default: |
| 108 | g_assert_not_reached(); |
| 109 | } |
| 110 | |
| 111 | if (is_snan(ret->cls)) { |
| 112 | return partsN(silence_nan)(ret, s); |
| 113 | } |
| 114 | return *ret; |
| 115 | } |
| 116 | |
| 117 | static FloatPartsN partsN(pick_nan_muladd)(const FloatPartsN *a, |
| 118 | const FloatPartsN *b, |
| 119 | const FloatPartsN *c, |
| 120 | float_status *s, |
| 121 | int ab_mask, int abc_mask) |
| 122 | { |
| 123 | bool infzero = (ab_mask == float_cmask_infzero); |
| 124 | bool have_snan = (abc_mask & float_cmask_snan); |
| 125 | FloatInfZeroNaNRule izn_rule = get_float_infzeronan_rule(s); |
| 126 | const FloatPartsN *ret; |
| 127 | |
| 128 | if (unlikely(have_snan)) { |
| 129 | float_raise(float_flag_invalid | float_flag_invalid_snan, s); |
| 130 | } |
| 131 | |
| 132 | if (infzero && !(izn_rule & float_infzeronan_suppress_invalid)) { |
| 133 | /* This is (0 * inf) + NaN or (inf * 0) + NaN */ |
| 134 | float_raise(float_flag_invalid | float_flag_invalid_imz, s); |
| 135 | } |
| 136 | |
| 137 | if (get_default_nan_mode(s)) { |
| 138 | /* |
| 139 | * We guarantee not to require the target to tell us how to |
| 140 | * pick a NaN if we're always returning the default NaN. |
| 141 | * But if we're not in default-NaN mode then the target must |
| 142 | * specify. |
| 143 | */ |
| 144 | goto default_nan; |
| 145 | } else if (infzero) { |
| 146 | /* |
| 147 | * Inf * 0 + NaN -- some implementations return the |
| 148 | * default NaN here, and some return the input NaN. |
| 149 | */ |
| 150 | switch (izn_rule & ~float_infzeronan_suppress_invalid) { |
| 151 | case float_infzeronan_dnan_never: |
| 152 | break; |
| 153 | case float_infzeronan_dnan_always: |
| 154 | goto default_nan; |
| 155 | case float_infzeronan_dnan_if_qnan: |
| 156 | if (is_qnan(c->cls)) { |
| 157 | goto default_nan; |
| 158 | } |
| 159 | break; |
| 160 | default: |
| 161 | g_assert_not_reached(); |
| 162 | } |
| 163 | ret = c; |
| 164 | } else { |
| 165 | const FloatPartsN *val[R_3NAN_1ST_MASK + 1] = { a, b, c }; |
| 166 | Float3NaNPropRule rule = get_float_3nan_prop_rule(s); |
| 167 | |
| 168 | assert(rule != float_3nan_prop_none); |
| 169 | if (have_snan && (rule & R_3NAN_SNAN_MASK)) { |
| 170 | /* We have at least one SNaN input and should prefer it */ |
| 171 | do { |
| 172 | ret = val[rule & R_3NAN_1ST_MASK]; |
| 173 | rule >>= R_3NAN_1ST_LENGTH; |
| 174 | } while (!is_snan(ret->cls)); |
| 175 | } else { |
| 176 | do { |
| 177 | ret = val[rule & R_3NAN_1ST_MASK]; |
| 178 | rule >>= R_3NAN_1ST_LENGTH; |
| 179 | } while (!is_nan(ret->cls)); |
| 180 | } |
| 181 | } |
| 182 | |
| 183 | if (is_snan(ret->cls)) { |
| 184 | return partsN(silence_nan)(ret, s); |
| 185 | } |
| 186 | return *ret; |
| 187 | |
| 188 | default_nan: |
| 189 | return partsN(default_nan)(s); |
| 190 | } |
| 191 | |
| 192 | /* |
| 193 | * Canonicalize the FloatParts structure. Determine the class, |
| 194 | * unbias the exponent, and normalize the fraction. |
| 195 | */ |
| 196 | static void partsN(canonicalize)(FloatPartsN *p, float_status *status, |
| 197 | const FloatFmt *fmt) |
| 198 | { |
| 199 | /* |
| 200 | * It's target-dependent how to handle the case of exponent 0 |
| 201 | * and Integer bit set. Intel calls these "pseudodenormals", |
| 202 | * and treats them as if the integer bit was 0, and never |
| 203 | * produces them on output. This is the default behaviour for QEMU. |
| 204 | * For m68k, the integer bit is considered validly part of the |
| 205 | * input value when the exponent is 0, and may be 0 or 1, |
| 206 | * giving extra range. They may also be generated as outputs. |
| 207 | * (The m68k manual actually calls these values part of the |
| 208 | * normalized number range, not the denormalized number range, |
| 209 | * but that distinction is not important for us, because |
| 210 | * m68k doesn't care about the input_denormal_used status flag.) |
| 211 | * floatx80_pseudo_denormal_valid selects the m68k behaviour, |
| 212 | * which changes both how we canonicalize such a value and |
| 213 | * how we uncanonicalize results. |
| 214 | */ |
| 215 | bool has_pseudo_denormals = fmt->has_explicit_bit && |
| 216 | (get_floatx80_behaviour(status) & floatx80_pseudo_denormal_valid); |
| 217 | |
| 218 | if (unlikely(p->exp == 0)) { |
| 219 | if (likely(fracN(eqz)(p))) { |
| 220 | p->cls = float_class_zero; |
| 221 | } else if (get_flush_inputs_to_zero(status)) { |
| 222 | float_raise(float_flag_input_denormal_flushed, status); |
| 223 | p->cls = float_class_zero; |
| 224 | fracN(clear)(p); |
| 225 | } else { |
| 226 | int shift = fracN(normalize)(p); |
| 227 | p->cls = float_class_denormal; |
| 228 | p->exp = fmt->frac_shift - fmt->exp_bias |
| 229 | - shift + !has_pseudo_denormals; |
| 230 | } |
| 231 | return; |
| 232 | } |
| 233 | if (unlikely(p->exp == fmt->exp_max)) { |
| 234 | switch (fmt->exp_max_kind) { |
| 235 | case float_expmax_ieee: |
| 236 | if (likely(fracN(eqz)(p))) { |
| 237 | p->cls = float_class_inf; |
| 238 | } else { |
| 239 | fracN(shl)(p, fmt->frac_shift); |
| 240 | p->cls = (parts_is_snan_frac(p->frac_hi, status) |
| 241 | ? float_class_snan : float_class_qnan); |
| 242 | } |
| 243 | return; |
| 244 | case float_expmax_normal: |
| 245 | break; |
| 246 | case float_expmax_e4m3: |
| 247 | if (p->frac_hi == 0b111) { |
| 248 | fracN(shl)(p, fmt->frac_shift); |
| 249 | p->cls = (get_float_e4m3_nan_is_snan(status) |
| 250 | ? float_class_snan : float_class_qnan); |
| 251 | return; |
| 252 | } |
| 253 | /* otherwise normal */ |
| 254 | break; |
| 255 | default: |
| 256 | g_assert_not_reached(); |
| 257 | } |
| 258 | } |
| 259 | |
| 260 | p->cls = float_class_normal; |
| 261 | p->exp -= fmt->exp_bias; |
| 262 | fracN(shl)(p, fmt->frac_shift); |
| 263 | p->frac_hi |= DECOMPOSED_IMPLICIT_BIT; |
| 264 | } |
| 265 | |
| 266 | /* |
| 267 | * Round and uncanonicalize a floating-point number by parts. There |
| 268 | * are FRAC_SHIFT bits that may require rounding at the bottom of the |
| 269 | * fraction; these bits will be removed. The exponent will be biased |
| 270 | * by EXP_BIAS and must be bounded by [EXP_MAX-1, 0]. |
| 271 | * |
| 272 | * The saturate parameter controls saturation behavior for formats that |
| 273 | * support it -- when true, overflow produces max normal instead of infinity. |
| 274 | */ |
| 275 | |
| 276 | /* Helper for uncanon_normal and uncanon, for FP8 E4M3. */ |
| 277 | static void partsN(uncanon_e4m3_overflow)(FloatPartsN *p, float_status *s, |
| 278 | const FloatFmt *fmt, bool saturate) |
| 279 | { |
| 280 | assert(N == 64); |
| 281 | p->exp = fmt->exp_max; |
| 282 | if (saturate) { |
| 283 | p->frac_hi = E4M3_NORMAL_FRAC_MAX; |
| 284 | } else { |
| 285 | /* |
| 286 | * The class isn't actually used after this point in uncanon, |
| 287 | * but for clarity while debugging, don't leave it set to normal. |
| 288 | */ |
| 289 | p->cls = float_class_qnan; |
| 290 | p->frac_hi = E4M3_NAN_FRAC; |
| 291 | } |
| 292 | } |
| 293 | |
| 294 | static void partsN(uncanon_normal)(FloatPartsN *p, float_status *s, |
| 295 | const FloatFmt *fmt, bool saturate) |
| 296 | { |
| 297 | const int exp_max = fmt->exp_max; |
| 298 | const int frac_shift = fmt->frac_shift; |
| 299 | const uint64_t round_mask = fmt->round_mask; |
| 300 | const uint64_t frac_lsb = round_mask + 1; |
| 301 | const uint64_t frac_lsbm1 = round_mask ^ (round_mask >> 1); |
| 302 | const uint64_t roundeven_mask = round_mask | frac_lsb; |
| 303 | uint64_t inc; |
| 304 | bool overflow_norm = saturate; |
| 305 | int exp; |
| 306 | FloatExceptionFlags flags = 0; |
| 307 | |
| 308 | switch (get_float_rounding_mode(s)) { |
| 309 | case float_round_nearest_even_max: |
| 310 | overflow_norm = true; |
| 311 | /* fall through */ |
| 312 | case float_round_nearest_even: |
| 313 | if (N > 64 && frac_lsb == 0) { |
| 314 | inc = ((p->frac_hi & 1) || (p->frac_lo & round_mask) != frac_lsbm1 |
| 315 | ? frac_lsbm1 : 0); |
| 316 | } else { |
| 317 | inc = ((p->frac_lo & roundeven_mask) != frac_lsbm1 |
| 318 | ? frac_lsbm1 : 0); |
| 319 | } |
| 320 | break; |
| 321 | case float_round_ties_away: |
| 322 | inc = frac_lsbm1; |
| 323 | break; |
| 324 | case float_round_to_zero: |
| 325 | overflow_norm = true; |
| 326 | inc = 0; |
| 327 | break; |
| 328 | case float_round_up: |
| 329 | inc = p->sign ? 0 : round_mask; |
| 330 | overflow_norm |= p->sign; |
| 331 | break; |
| 332 | case float_round_down: |
| 333 | inc = p->sign ? round_mask : 0; |
| 334 | overflow_norm |= !p->sign; |
| 335 | break; |
| 336 | case float_round_to_odd: |
| 337 | overflow_norm = true; |
| 338 | /* fall through */ |
| 339 | case float_round_to_odd_inf: |
| 340 | if (N > 64 && frac_lsb == 0) { |
| 341 | inc = p->frac_hi & 1 ? 0 : round_mask; |
| 342 | } else { |
| 343 | inc = p->frac_lo & frac_lsb ? 0 : round_mask; |
| 344 | } |
| 345 | break; |
| 346 | default: |
| 347 | g_assert_not_reached(); |
| 348 | } |
| 349 | |
| 350 | exp = p->exp + fmt->exp_bias; |
| 351 | if (likely(exp > 0)) { |
| 352 | if (p->frac_lo & round_mask) { |
| 353 | flags |= float_flag_inexact; |
| 354 | if (fracN(addi)(p, p, inc)) { |
| 355 | fracN(shr)(p, 1); |
| 356 | p->frac_hi |= DECOMPOSED_IMPLICIT_BIT; |
| 357 | exp++; |
| 358 | } |
| 359 | p->frac_lo &= ~round_mask; |
| 360 | } |
| 361 | |
| 362 | if (unlikely(exp >= exp_max)) { |
| 363 | switch (fmt->exp_max_kind) { |
| 364 | case float_expmax_ieee: |
| 365 | flags |= float_flag_overflow; |
| 366 | if (get_float_rebias_overflow(s)) { |
| 367 | exp -= fmt->exp_re_bias; |
| 368 | } else if (overflow_norm) { |
| 369 | flags |= float_flag_inexact; |
| 370 | exp = exp_max - 1; |
| 371 | fracN(allones)(p); |
| 372 | p->frac_lo &= ~round_mask; |
| 373 | } else { |
| 374 | flags |= float_flag_inexact; |
| 375 | p->cls = float_class_inf; |
| 376 | exp = exp_max; |
| 377 | fracN(clear)(p); |
| 378 | } |
| 379 | break; |
| 380 | |
| 381 | case float_expmax_normal: |
| 382 | if (unlikely(exp > exp_max)) { |
| 383 | /* Overflow. Return the maximum normal. */ |
| 384 | flags = (fmt->overflow_raises_invalid |
| 385 | ? float_flag_invalid |
| 386 | : float_flag_overflow | float_flag_inexact); |
| 387 | exp = exp_max; |
| 388 | fracN(allones)(p); |
| 389 | p->frac_lo &= ~round_mask; |
| 390 | } |
| 391 | break; |
| 392 | |
| 393 | case float_expmax_e4m3: |
| 394 | if (exp > exp_max || p->frac_hi > E4M3_NORMAL_FRAC_MAX) { |
| 395 | partsN(uncanon_e4m3_overflow)(p, s, fmt, overflow_norm); |
| 396 | exp = p->exp; |
| 397 | flags |= (float_flag_overflow | float_flag_inexact); |
| 398 | } |
| 399 | break; |
| 400 | |
| 401 | default: |
| 402 | g_assert_not_reached(); |
| 403 | } |
| 404 | } |
| 405 | fracN(shr)(p, frac_shift); |
| 406 | } else if (unlikely(get_float_rebias_underflow(s))) { |
| 407 | flags |= float_flag_underflow; |
| 408 | exp += fmt->exp_re_bias; |
| 409 | if (p->frac_lo & round_mask) { |
| 410 | flags |= float_flag_inexact; |
| 411 | if (fracN(addi)(p, p, inc)) { |
| 412 | fracN(shr)(p, 1); |
| 413 | p->frac_hi |= DECOMPOSED_IMPLICIT_BIT; |
| 414 | exp++; |
| 415 | } |
| 416 | p->frac_lo &= ~round_mask; |
| 417 | } |
| 418 | fracN(shr)(p, frac_shift); |
| 419 | } else if (get_flush_to_zero(s) && get_ftz_before_rounding(s)) { |
| 420 | flags |= float_flag_output_denormal_flushed; |
| 421 | p->cls = float_class_zero; |
| 422 | exp = 0; |
| 423 | fracN(clear)(p); |
| 424 | } else { |
| 425 | bool is_tiny = get_tininess_before_rounding(s) || exp < 0; |
| 426 | bool has_pseudo_denormals = fmt->has_explicit_bit && |
| 427 | (get_floatx80_behaviour(s) & floatx80_pseudo_denormal_valid); |
| 428 | |
| 429 | if (!is_tiny) { |
| 430 | FloatPartsN discard; |
| 431 | is_tiny = !fracN(addi)(&discard, p, inc); |
| 432 | } |
| 433 | |
| 434 | fracN(shrjam)(p, !has_pseudo_denormals - exp); |
| 435 | |
| 436 | if (p->frac_lo & round_mask) { |
| 437 | /* Need to recompute round-to-even/round-to-odd. */ |
| 438 | switch (get_float_rounding_mode(s)) { |
| 439 | case float_round_nearest_even: |
| 440 | case float_round_nearest_even_max: |
| 441 | if (N > 64 && frac_lsb == 0) { |
| 442 | inc = ((p->frac_hi & 1) || |
| 443 | (p->frac_lo & round_mask) != frac_lsbm1 |
| 444 | ? frac_lsbm1 : 0); |
| 445 | } else { |
| 446 | inc = ((p->frac_lo & roundeven_mask) != frac_lsbm1 |
| 447 | ? frac_lsbm1 : 0); |
| 448 | } |
| 449 | break; |
| 450 | case float_round_to_odd: |
| 451 | case float_round_to_odd_inf: |
| 452 | if (N > 64 && frac_lsb == 0) { |
| 453 | inc = p->frac_hi & 1 ? 0 : round_mask; |
| 454 | } else { |
| 455 | inc = p->frac_lo & frac_lsb ? 0 : round_mask; |
| 456 | } |
| 457 | break; |
| 458 | default: |
| 459 | break; |
| 460 | } |
| 461 | flags |= float_flag_inexact; |
| 462 | fracN(addi)(p, p, inc); |
| 463 | p->frac_lo &= ~round_mask; |
| 464 | } |
| 465 | |
| 466 | exp = (p->frac_hi & DECOMPOSED_IMPLICIT_BIT) && !has_pseudo_denormals; |
| 467 | fracN(shr)(p, frac_shift); |
| 468 | |
| 469 | if (is_tiny) { |
| 470 | if (get_flush_to_zero(s)) { |
| 471 | assert(!get_ftz_before_rounding(s)); |
| 472 | flags |= float_flag_output_denormal_flushed; |
| 473 | p->cls = float_class_zero; |
| 474 | exp = 0; |
| 475 | fracN(clear)(p); |
| 476 | } else if (flags & float_flag_inexact) { |
| 477 | flags |= float_flag_underflow; |
| 478 | } |
| 479 | if (exp == 0 && fracN(eqz)(p)) { |
| 480 | p->cls = float_class_zero; |
| 481 | } |
| 482 | } |
| 483 | } |
| 484 | p->exp = exp; |
| 485 | float_raise(flags, s); |
| 486 | } |
| 487 | |
| 488 | static void partsN(uncanon)(FloatPartsN *p, float_status *s, |
| 489 | const FloatFmt *fmt, bool saturate) |
| 490 | { |
| 491 | if (likely(is_anynorm(p->cls))) { |
| 492 | partsN(uncanon_normal)(p, s, fmt, saturate); |
| 493 | } else { |
| 494 | switch (p->cls) { |
| 495 | case float_class_zero: |
| 496 | p->exp = 0; |
| 497 | fracN(clear)(p); |
| 498 | return; |
| 499 | case float_class_inf: |
| 500 | switch (fmt->exp_max_kind) { |
| 501 | case float_expmax_ieee: |
| 502 | p->exp = fmt->exp_max; |
| 503 | fracN(clear)(p); |
| 504 | break; |
| 505 | case float_expmax_e4m3: |
| 506 | partsN(uncanon_e4m3_overflow)(p, s, fmt, saturate); |
| 507 | fracN(shr)(p, fmt->frac_shift); |
| 508 | break; |
| 509 | case float_expmax_normal: |
| 510 | default: |
| 511 | g_assert_not_reached(); |
| 512 | } |
| 513 | return; |
| 514 | case float_class_qnan: |
| 515 | case float_class_snan: |
| 516 | p->exp = fmt->exp_max; |
| 517 | switch (fmt->exp_max_kind) { |
| 518 | case float_expmax_e4m3: |
| 519 | /* |
| 520 | * There is only one NaN encoding for E4M3, and with a |
| 521 | * conversion from another format, the input NaN fraction |
| 522 | * may not apply. |
| 523 | */ |
| 524 | assert(N == 64); |
| 525 | p->frac_hi = E4M3_NAN_FRAC; |
| 526 | /* fall through */ |
| 527 | case float_expmax_ieee: |
| 528 | fracN(shr)(p, fmt->frac_shift); |
| 529 | break; |
| 530 | case float_expmax_normal: |
| 531 | default: |
| 532 | g_assert_not_reached(); |
| 533 | } |
| 534 | return; |
| 535 | default: |
| 536 | break; |
| 537 | } |
| 538 | g_assert_not_reached(); |
| 539 | } |
| 540 | } |
| 541 | |
| 542 | /* |
| 543 | * Returns the result of adding or subtracting the values of the |
| 544 | * floating-point values `a' and `b'. The operation is performed |
| 545 | * according to the IEC/IEEE Standard for Binary Floating-Point |
| 546 | * Arithmetic. |
| 547 | */ |
| 548 | FloatPartsN partsN(addsub)(const FloatPartsN *a_orig, |
| 549 | const FloatPartsN *b_orig, |
| 550 | float_status *s, bool subtract) |
| 551 | { |
| 552 | int ab_mask = float_cmask(a_orig->cls) | float_cmask(b_orig->cls); |
| 553 | |
| 554 | if (unlikely(ab_mask & float_cmask_anynan)) { |
| 555 | return partsN(pick_nan)(a_orig, b_orig, s); |
| 556 | } |
| 557 | |
| 558 | /* |
| 559 | * For addition and subtraction, we will consume an |
| 560 | * input denormal unless the other input is a NaN. |
| 561 | */ |
| 562 | record_denormals_used(ab_mask, s); |
| 563 | |
| 564 | FloatPartsN a = *a_orig; |
| 565 | FloatPartsN b = *b_orig; |
| 566 | |
| 567 | b.sign ^= subtract; |
| 568 | |
| 569 | if (a.sign != b.sign) { |
| 570 | /* Subtraction */ |
| 571 | if (likely(cmask_is_only_normals(ab_mask))) { |
| 572 | if (partsN(sub_normal)(&a, &b)) { |
| 573 | return a; |
| 574 | } |
| 575 | /* Subtract was exact, fall through to set sign. */ |
| 576 | ab_mask = float_cmask_zero; |
| 577 | } |
| 578 | |
| 579 | if (ab_mask == float_cmask_zero) { |
| 580 | /* 0 - 0 */ |
| 581 | a.sign = get_float_rounding_mode(s) == float_round_down; |
| 582 | return a; |
| 583 | } |
| 584 | |
| 585 | if (ab_mask & float_cmask_inf) { |
| 586 | if (a.cls != float_class_inf) { |
| 587 | /* N - Inf */ |
| 588 | return b; |
| 589 | } |
| 590 | if (b.cls != float_class_inf) { |
| 591 | /* Inf - N */ |
| 592 | return a; |
| 593 | } |
| 594 | /* Inf - Inf */ |
| 595 | float_raise(float_flag_invalid | float_flag_invalid_isi, s); |
| 596 | return partsN(default_nan)(s); |
| 597 | } |
| 598 | } else { |
| 599 | /* Addition */ |
| 600 | if (likely(cmask_is_only_normals(ab_mask))) { |
| 601 | partsN(add_normal)(&a, &b); |
| 602 | return a; |
| 603 | } |
| 604 | |
| 605 | if (ab_mask == float_cmask_zero) { |
| 606 | /* 0 + 0 */ |
| 607 | return a; |
| 608 | } |
| 609 | |
| 610 | if (ab_mask & float_cmask_inf) { |
| 611 | /* N + Inf or Inf + N */ |
| 612 | a.cls = float_class_inf; |
| 613 | return a; |
| 614 | } |
| 615 | } |
| 616 | |
| 617 | /* 0 +/- N or N +/- 0 */ |
| 618 | assert((ab_mask & float_cmask_zero) && (ab_mask & float_cmask_anynorm)); |
| 619 | return b.cls == float_class_zero ? a : b; |
| 620 | } |
| 621 | |
| 622 | /* |
| 623 | * Returns the result of multiplying the floating-point values `a' and |
| 624 | * `b'. The operation is performed according to the IEC/IEEE Standard |
| 625 | * for Binary Floating-Point Arithmetic. |
| 626 | */ |
| 627 | FloatPartsN partsN(mul)(const FloatPartsN *a, const FloatPartsN *b, |
| 628 | float_status *s) |
| 629 | { |
| 630 | int ab_mask = float_cmask(a->cls) | float_cmask(b->cls); |
| 631 | bool sign = a->sign ^ b->sign; |
| 632 | |
| 633 | if (likely(cmask_is_only_normals(ab_mask))) { |
| 634 | FloatPartsW tmp; |
| 635 | FloatPartsN r = { |
| 636 | .cls = float_class_normal, |
| 637 | .sign = sign, |
| 638 | .exp = a->exp + b->exp + 1, |
| 639 | }; |
| 640 | |
| 641 | record_denormals_used(ab_mask, s); |
| 642 | |
| 643 | fracN(mulw)(&tmp, a, b); |
| 644 | fracN(truncjam)(&r, &tmp); |
| 645 | |
| 646 | if (!(r.frac_hi & DECOMPOSED_IMPLICIT_BIT)) { |
| 647 | fracN(add)(&r, &r, &r); |
| 648 | r.exp -= 1; |
| 649 | } |
| 650 | |
| 651 | return r; |
| 652 | } |
| 653 | |
| 654 | /* Inf * Zero == NaN */ |
| 655 | if (unlikely(ab_mask == float_cmask_infzero)) { |
| 656 | float_raise(float_flag_invalid | float_flag_invalid_imz, s); |
| 657 | return partsN(default_nan)(s); |
| 658 | } |
| 659 | |
| 660 | if (unlikely(ab_mask & float_cmask_anynan)) { |
| 661 | return partsN(pick_nan)(a, b, s); |
| 662 | } |
| 663 | |
| 664 | /* Multiply by 0 or Inf */ |
| 665 | record_denormals_used(ab_mask, s); |
| 666 | |
| 667 | if (ab_mask & float_cmask_inf) { |
| 668 | return (FloatPartsN){ .cls = float_class_inf, .sign = sign }; |
| 669 | } |
| 670 | |
| 671 | g_assert(ab_mask & float_cmask_zero); |
| 672 | return (FloatPartsN){ .cls = float_class_zero, .sign = sign }; |
| 673 | } |
| 674 | |
| 675 | /* |
| 676 | * Returns the result of multiplying the floating-point values `a' and |
| 677 | * `b' then adding 'c', with no intermediate rounding step after the |
| 678 | * multiplication. The operation is performed according to the |
| 679 | * IEC/IEEE Standard for Binary Floating-Point Arithmetic 754-2008. |
| 680 | * The flags argument allows the caller to select negation of the addend |
| 681 | * or the intermediate product. (The difference between this and having |
| 682 | * the caller do a separate negation is that negating externally will |
| 683 | * flip the sign bit on NaNs.) Note that float_muladd_negate_result |
| 684 | * is not applied here, and should be handled separately after rounding |
| 685 | * chooses the final sign of 0.0. |
| 686 | * |
| 687 | * Requires A and C extracted into a double-sized structure to provide the |
| 688 | * extra space for the widening multiply. |
| 689 | */ |
| 690 | FloatPartsN partsN(muladd)(const FloatPartsN *a, const FloatPartsN *b, |
| 691 | const FloatPartsN *c, int flags, float_status *s) |
| 692 | { |
| 693 | int ab_mask = float_cmask(a->cls) | float_cmask(b->cls); |
| 694 | int c_mask = float_cmask(c->cls); |
| 695 | int abc_mask = ab_mask | c_mask; |
| 696 | bool c_sign = c->sign ^ !!(flags & float_muladd_negate_c); |
| 697 | bool p_sign = a->sign ^ b->sign ^ !!(flags & float_muladd_negate_product); |
| 698 | |
| 699 | /* |
| 700 | * The "likely" case is A and B normal, so that the product is normal, |
| 701 | * and C normal or zero so that the result is normal. |
| 702 | */ |
| 703 | int likely_mask = ab_mask | (c_mask & ~float_cmask_zero); |
| 704 | if (likely(cmask_is_only_normals(likely_mask))) { |
| 705 | record_denormals_used(abc_mask, s); |
| 706 | |
| 707 | /* Perform the multiplication step. */ |
| 708 | FloatPartsW p_widen = { .sign = p_sign, .exp = a->exp + b->exp + 1 }; |
| 709 | fracN(mulw)(&p_widen, a, b); |
| 710 | if (!(p_widen.frac_hi & DECOMPOSED_IMPLICIT_BIT)) { |
| 711 | fracW(add)(&p_widen, &p_widen, &p_widen); |
| 712 | p_widen.exp -= 1; |
| 713 | } |
| 714 | |
| 715 | /* Perform the addition step. */ |
| 716 | if (!(c_mask & float_cmask_zero)) { |
| 717 | /* Zero-extend C to less significant bits. */ |
| 718 | FloatPartsW c_widen = { .sign = c_sign, .exp = c->exp }; |
| 719 | fracN(widen)(&c_widen, c); |
| 720 | |
| 721 | if (p_sign == c_sign) { |
| 722 | partsW(add_normal)(&p_widen, &c_widen); |
| 723 | } else if (!partsW(sub_normal)(&p_widen, &c_widen)) { |
| 724 | goto return_sub_zero; |
| 725 | } |
| 726 | } |
| 727 | |
| 728 | /* Narrow with sticky bit, for proper rounding later. */ |
| 729 | FloatPartsN r = { |
| 730 | .sign = p_widen.sign, |
| 731 | .exp = p_widen.exp, |
| 732 | .cls = float_class_normal, |
| 733 | }; |
| 734 | fracN(truncjam)(&r, &p_widen); |
| 735 | return r; |
| 736 | } |
| 737 | |
| 738 | /* |
| 739 | * It is implementation-defined whether the cases of (0,inf,qnan) |
| 740 | * and (inf,0,qnan) raise InvalidOperation or not (and what QNaN |
| 741 | * they return if they do), so we have to hand this information |
| 742 | * off to the target-specific pick-a-NaN routine. |
| 743 | */ |
| 744 | if (unlikely(abc_mask & float_cmask_anynan)) { |
| 745 | return partsN(pick_nan_muladd)(a, b, c, s, ab_mask, abc_mask); |
| 746 | } |
| 747 | |
| 748 | if (unlikely(ab_mask == float_cmask_infzero)) { |
| 749 | /* Inf * Zero == NaN */ |
| 750 | float_raise(float_flag_invalid | float_flag_invalid_imz, s); |
| 751 | goto d_nan; |
| 752 | } |
| 753 | |
| 754 | if (unlikely(ab_mask & float_cmask_inf)) { |
| 755 | if ((c_mask & float_cmask_inf) && p_sign != c_sign) { |
| 756 | /* Inf - Inf == NaN */ |
| 757 | float_raise(float_flag_invalid | float_flag_invalid_isi, s); |
| 758 | goto d_nan; |
| 759 | } |
| 760 | /* Inf + C == Inf */ |
| 761 | record_denormals_used(abc_mask, s); |
| 762 | return (FloatPartsN){ .sign = p_sign, .cls = float_class_inf }; |
| 763 | } |
| 764 | record_denormals_used(abc_mask, s); |
| 765 | |
| 766 | /* Only remaining cases are zero product or inf addend. */ |
| 767 | assert((ab_mask & float_cmask_zero) | (c_mask & float_cmask_inf)); |
| 768 | |
| 769 | /* |
| 770 | * P + Inf == Inf, or |
| 771 | * 0 + C == C, |
| 772 | * except for 0 - 0, which needs special rounding, |
| 773 | * except for when we want to suppress this addition step. |
| 774 | */ |
| 775 | if (!(c_mask & float_cmask_zero) |
| 776 | || p_sign == c_sign |
| 777 | || (flags & float_muladd_suppress_add_product_zero)) { |
| 778 | FloatPartsN r = *c; |
| 779 | r.sign = c_sign; |
| 780 | return r; |
| 781 | } |
| 782 | |
| 783 | return_sub_zero: |
| 784 | /* 0 - 0 == -0 for round_down, +0 otherwise. */ |
| 785 | return (FloatPartsN){ |
| 786 | .sign = get_float_rounding_mode(s) == float_round_down, |
| 787 | .cls = float_class_zero |
| 788 | }; |
| 789 | |
| 790 | d_nan: |
| 791 | return partsN(default_nan)(s); |
| 792 | } |
| 793 | |
| 794 | /* |
| 795 | * Returns the result of dividing the floating-point value `a' by the |
| 796 | * corresponding value `b'. The operation is performed according to |
| 797 | * the IEC/IEEE Standard for Binary Floating-Point Arithmetic. |
| 798 | */ |
| 799 | FloatPartsN partsN(div)(const FloatPartsN *a, const FloatPartsN *b, |
| 800 | float_status *s) |
| 801 | { |
| 802 | int ab_mask = float_cmask(a->cls) | float_cmask(b->cls); |
| 803 | FloatPartsN r = *a; |
| 804 | |
| 805 | r.sign ^= b->sign; |
| 806 | r.exp -= b->exp; |
| 807 | |
| 808 | if (likely(cmask_is_only_normals(ab_mask))) { |
| 809 | record_denormals_used(ab_mask, s); |
| 810 | r.exp -= fracN(div)(&r, b); |
| 811 | return r; |
| 812 | } |
| 813 | |
| 814 | /* 0/0 or Inf/Inf => NaN */ |
| 815 | if (unlikely(ab_mask == float_cmask_zero)) { |
| 816 | float_raise(float_flag_invalid | float_flag_invalid_zdz, s); |
| 817 | return partsN(default_nan)(s); |
| 818 | } |
| 819 | if (unlikely(ab_mask == float_cmask_inf)) { |
| 820 | float_raise(float_flag_invalid | float_flag_invalid_idi, s); |
| 821 | return partsN(default_nan)(s); |
| 822 | } |
| 823 | |
| 824 | /* All the NaN cases */ |
| 825 | if (unlikely(ab_mask & float_cmask_anynan)) { |
| 826 | return partsN(pick_nan)(a, b, s); |
| 827 | } |
| 828 | |
| 829 | if (b->cls != float_class_zero) { |
| 830 | record_denormals_used(ab_mask, s); |
| 831 | } |
| 832 | |
| 833 | /* Inf / X */ |
| 834 | if (r.cls == float_class_inf) { |
| 835 | return r; |
| 836 | } |
| 837 | |
| 838 | /* 0 / X */ |
| 839 | if (r.cls == float_class_zero) { |
| 840 | return r; |
| 841 | } |
| 842 | |
| 843 | /* X / Inf */ |
| 844 | if (b->cls == float_class_inf) { |
| 845 | r.cls = float_class_zero; |
| 846 | return r; |
| 847 | } |
| 848 | |
| 849 | /* X / 0 => Inf */ |
| 850 | assert(b->cls == float_class_zero); |
| 851 | float_raise(float_flag_divbyzero, s); |
| 852 | r.cls = float_class_inf; |
| 853 | return r; |
| 854 | } |
| 855 | |
| 856 | /* |
| 857 | * Floating point remainder, per IEC/IEEE, or modulus. |
| 858 | */ |
| 859 | static FloatPartsN *partsN(modrem)(FloatPartsN *a, FloatPartsN *b, |
| 860 | uint64_t *mod_quot, float_status *s) |
| 861 | { |
| 862 | int ab_mask = float_cmask(a->cls) | float_cmask(b->cls); |
| 863 | |
| 864 | if (likely(cmask_is_only_normals(ab_mask))) { |
| 865 | record_denormals_used(ab_mask, s); |
| 866 | fracN(modrem)(a, b, mod_quot); |
| 867 | return a; |
| 868 | } |
| 869 | |
| 870 | if (mod_quot) { |
| 871 | *mod_quot = 0; |
| 872 | } |
| 873 | |
| 874 | /* All the NaN cases */ |
| 875 | if (unlikely(ab_mask & float_cmask_anynan)) { |
| 876 | *a = partsN(pick_nan)(a, b, s); |
| 877 | return a; |
| 878 | } |
| 879 | |
| 880 | /* Inf % N; N % 0 */ |
| 881 | if (a->cls == float_class_inf || b->cls == float_class_zero) { |
| 882 | float_raise(float_flag_invalid, s); |
| 883 | *a = partsN(default_nan)(s); |
| 884 | return a; |
| 885 | } |
| 886 | |
| 887 | record_denormals_used(ab_mask, s); |
| 888 | |
| 889 | /* N % Inf; 0 % N */ |
| 890 | g_assert(b->cls == float_class_inf || a->cls == float_class_zero); |
| 891 | return a; |
| 892 | } |
| 893 | |
| 894 | /* |
| 895 | * Square Root |
| 896 | * |
| 897 | * The base algorithm is lifted from |
| 898 | * https://git.musl-libc.org/cgit/musl/tree/src/math/sqrtf.c |
| 899 | * https://git.musl-libc.org/cgit/musl/tree/src/math/sqrt.c |
| 900 | * https://git.musl-libc.org/cgit/musl/tree/src/math/sqrtl.c |
| 901 | * and is thus MIT licenced. |
| 902 | */ |
| 903 | static void partsN(sqrt)(FloatPartsN *a, float_status *status, |
| 904 | const FloatFmt *fmt) |
| 905 | { |
| 906 | const uint32_t three32 = 3u << 30; |
| 907 | const uint64_t three64 = 3ull << 62; |
| 908 | uint32_t d32, m32, r32, s32, u32; /* 32-bit computation */ |
| 909 | uint64_t d64, m64, r64, s64, u64; /* 64-bit computation */ |
| 910 | uint64_t dh, dl, rh, rl, sh, sl, uh, ul; /* 128-bit computation */ |
| 911 | uint64_t d0h, d0l, d1h, d1l, d2h, d2l; |
| 912 | uint64_t discard; |
| 913 | bool exp_odd; |
| 914 | size_t index; |
| 915 | |
| 916 | if (unlikely(a->cls != float_class_normal)) { |
| 917 | switch (a->cls) { |
| 918 | case float_class_denormal: |
| 919 | if (!a->sign) { |
| 920 | /* -ve denormal will be InvalidOperation */ |
| 921 | float_raise(float_flag_input_denormal_used, status); |
| 922 | } |
| 923 | break; |
| 924 | case float_class_snan: |
| 925 | case float_class_qnan: |
| 926 | *a = partsN(return_nan)(a, status); |
| 927 | return; |
| 928 | case float_class_zero: |
| 929 | return; |
| 930 | case float_class_inf: |
| 931 | if (unlikely(a->sign)) { |
| 932 | goto d_nan; |
| 933 | } |
| 934 | return; |
| 935 | default: |
| 936 | g_assert_not_reached(); |
| 937 | } |
| 938 | } |
| 939 | |
| 940 | if (unlikely(a->sign)) { |
| 941 | goto d_nan; |
| 942 | } |
| 943 | |
| 944 | /* |
| 945 | * Argument reduction. |
| 946 | * x = 4^e frac; with integer e, and frac in [1, 4) |
| 947 | * m = frac fixed point at bit 62, since we're in base 4. |
| 948 | * If base-2 exponent is odd, exchange that for multiply by 2, |
| 949 | * which results in no shift. |
| 950 | */ |
| 951 | exp_odd = a->exp & 1; |
| 952 | index = extract64(a->frac_hi, 57, 6) | (!exp_odd << 6); |
| 953 | if (!exp_odd) { |
| 954 | fracN(shr)(a, 1); |
| 955 | } |
| 956 | |
| 957 | /* |
| 958 | * Approximate r ~= 1/sqrt(m) and s ~= sqrt(m) when m in [1, 4). |
| 959 | * |
| 960 | * Initial estimate: |
| 961 | * 7-bit lookup table (1-bit exponent and 6-bit significand). |
| 962 | * |
| 963 | * The relative error (e = r0*sqrt(m)-1) of a linear estimate |
| 964 | * (r0 = a*m + b) is |e| < 0.085955 ~ 0x1.6p-4 at best; |
| 965 | * a table lookup is faster and needs one less iteration. |
| 966 | * The 7-bit table gives |e| < 0x1.fdp-9. |
| 967 | * |
| 968 | * A Newton-Raphson iteration for r is |
| 969 | * s = m*r |
| 970 | * d = s*r |
| 971 | * u = 3 - d |
| 972 | * r = r*u/2 |
| 973 | * |
| 974 | * Fixed point representations: |
| 975 | * m, s, d, u, three are all 2.30; r is 0.32 |
| 976 | */ |
| 977 | m64 = a->frac_hi; |
| 978 | m32 = m64 >> 32; |
| 979 | |
| 980 | r32 = rsqrt_tab[index] << 16; |
| 981 | /* |r*sqrt(m) - 1| < 0x1.FDp-9 */ |
| 982 | |
| 983 | s32 = ((uint64_t)m32 * r32) >> 32; |
| 984 | d32 = ((uint64_t)s32 * r32) >> 32; |
| 985 | u32 = three32 - d32; |
| 986 | |
| 987 | if (N == 64) { |
| 988 | /* float64 or smaller */ |
| 989 | |
| 990 | r32 = ((uint64_t)r32 * u32) >> 31; |
| 991 | /* |r*sqrt(m) - 1| < 0x1.7Bp-16 */ |
| 992 | |
| 993 | s32 = ((uint64_t)m32 * r32) >> 32; |
| 994 | d32 = ((uint64_t)s32 * r32) >> 32; |
| 995 | u32 = three32 - d32; |
| 996 | |
| 997 | if (fmt->frac_size <= 23) { |
| 998 | /* float32 or smaller */ |
| 999 | |
| 1000 | s32 = ((uint64_t)s32 * u32) >> 32; /* 3.29 */ |
| 1001 | s32 = (s32 - 1) >> 6; /* 9.23 */ |
| 1002 | /* s < sqrt(m) < s + 0x1.08p-23 */ |
| 1003 | |
| 1004 | /* compute nearest rounded result to 2.23 bits */ |
| 1005 | uint32_t d0 = (m32 << 16) - s32 * s32; |
| 1006 | uint32_t d1 = s32 - d0; |
| 1007 | uint32_t d2 = d1 + s32 + 1; |
| 1008 | s32 += d1 >> 31; |
| 1009 | a->frac_hi = (uint64_t)s32 << (64 - 25); |
| 1010 | |
| 1011 | /* increment or decrement for inexact */ |
| 1012 | if (d2 != 0) { |
| 1013 | a->frac_hi += ((int32_t)(d1 ^ d2) < 0 ? -1 : 1); |
| 1014 | } |
| 1015 | goto done; |
| 1016 | } |
| 1017 | |
| 1018 | /* float64 */ |
| 1019 | |
| 1020 | r64 = (uint64_t)r32 * u32 * 2; |
| 1021 | /* |r*sqrt(m) - 1| < 0x1.37-p29; convert to 64-bit arithmetic */ |
| 1022 | mul64To128(m64, r64, &s64, &discard); |
| 1023 | mul64To128(s64, r64, &d64, &discard); |
| 1024 | u64 = three64 - d64; |
| 1025 | |
| 1026 | mul64To128(s64, u64, &s64, &discard); /* 3.61 */ |
| 1027 | s64 = (s64 - 2) >> 9; /* 12.52 */ |
| 1028 | |
| 1029 | /* Compute nearest rounded result */ |
| 1030 | uint64_t d0 = (m64 << 42) - s64 * s64; |
| 1031 | uint64_t d1 = s64 - d0; |
| 1032 | uint64_t d2 = d1 + s64 + 1; |
| 1033 | s64 += d1 >> 63; |
| 1034 | a->frac_hi = s64 << (64 - 54); |
| 1035 | |
| 1036 | /* increment or decrement for inexact */ |
| 1037 | if (d2 != 0) { |
| 1038 | a->frac_hi += ((int64_t)(d1 ^ d2) < 0 ? -1 : 1); |
| 1039 | } |
| 1040 | goto done; |
| 1041 | } |
| 1042 | |
| 1043 | r64 = (uint64_t)r32 * u32 * 2; |
| 1044 | /* |r*sqrt(m) - 1| < 0x1.7Bp-16; convert to 64-bit arithmetic */ |
| 1045 | |
| 1046 | mul64To128(m64, r64, &s64, &discard); |
| 1047 | mul64To128(s64, r64, &d64, &discard); |
| 1048 | u64 = three64 - d64; |
| 1049 | mul64To128(u64, r64, &r64, &discard); |
| 1050 | r64 <<= 1; |
| 1051 | /* |r*sqrt(m) - 1| < 0x1.a5p-31 */ |
| 1052 | |
| 1053 | mul64To128(m64, r64, &s64, &discard); |
| 1054 | mul64To128(s64, r64, &d64, &discard); |
| 1055 | u64 = three64 - d64; |
| 1056 | mul64To128(u64, r64, &rh, &rl); |
| 1057 | add128(rh, rl, rh, rl, &rh, &rl); |
| 1058 | /* |r*sqrt(m) - 1| < 0x1.c001p-59; change to 128-bit arithmetic */ |
| 1059 | |
| 1060 | mul128To256(a->frac_hi, a->frac_lo, rh, rl, &sh, &sl, &discard, &discard); |
| 1061 | mul128To256(sh, sl, rh, rl, &dh, &dl, &discard, &discard); |
| 1062 | sub128(three64, 0, dh, dl, &uh, &ul); |
| 1063 | mul128To256(uh, ul, sh, sl, &sh, &sl, &discard, &discard); /* 3.125 */ |
| 1064 | /* -0x1p-116 < s - sqrt(m) < 0x3.8001p-125 */ |
| 1065 | |
| 1066 | sub128(sh, sl, 0, 4, &sh, &sl); |
| 1067 | shift128Right(sh, sl, 13, &sh, &sl); /* 16.112 */ |
| 1068 | /* s < sqrt(m) < s + 1ulp */ |
| 1069 | |
| 1070 | /* Compute nearest rounded result */ |
| 1071 | mul64To128(sl, sl, &d0h, &d0l); |
| 1072 | d0h += 2 * sh * sl; |
| 1073 | sub128(a->frac_lo << 34, 0, d0h, d0l, &d0h, &d0l); |
| 1074 | sub128(sh, sl, d0h, d0l, &d1h, &d1l); |
| 1075 | add128(sh, sl, 0, 1, &d2h, &d2l); |
| 1076 | add128(d2h, d2l, d1h, d1l, &d2h, &d2l); |
| 1077 | add128(sh, sl, 0, d1h >> 63, &sh, &sl); |
| 1078 | shift128Left(sh, sl, 128 - 114, &sh, &sl); |
| 1079 | |
| 1080 | /* increment or decrement for inexact */ |
| 1081 | if (d2h | d2l) { |
| 1082 | if ((int64_t)(d1h ^ d2h) < 0) { |
| 1083 | sub128(sh, sl, 0, 1, &sh, &sl); |
| 1084 | } else { |
| 1085 | add128(sh, sl, 0, 1, &sh, &sl); |
| 1086 | } |
| 1087 | } |
| 1088 | a->frac_lo = sl; |
| 1089 | a->frac_hi = sh; |
| 1090 | |
| 1091 | done: |
| 1092 | /* Convert back from base 4 to base 2. */ |
| 1093 | a->exp >>= 1; |
| 1094 | if (!(a->frac_hi & DECOMPOSED_IMPLICIT_BIT)) { |
| 1095 | fracN(add)(a, a, a); |
| 1096 | } else { |
| 1097 | a->exp += 1; |
| 1098 | } |
| 1099 | return; |
| 1100 | |
| 1101 | d_nan: |
| 1102 | float_raise(float_flag_invalid | float_flag_invalid_sqrt, status); |
| 1103 | *a = partsN(default_nan)(status); |
| 1104 | } |
| 1105 | |
| 1106 | /* |
| 1107 | * Rounds the floating-point value `a' to an integer, and returns the |
| 1108 | * result as a floating-point value. The operation is performed |
| 1109 | * according to the IEC/IEEE Standard for Binary Floating-Point |
| 1110 | * Arithmetic. |
| 1111 | * |
| 1112 | * partsN(round_to_int_normal) is an internal helper function for |
| 1113 | * normal numbers only, returning true for inexact but not directly |
| 1114 | * raising float_flag_inexact. |
| 1115 | */ |
| 1116 | static bool partsN(round_to_int_normal)(FloatPartsN *a, FloatRoundMode rmode, |
| 1117 | int scale, int frac_size) |
| 1118 | { |
| 1119 | uint64_t frac_lsb, frac_lsbm1, rnd_even_mask, rnd_mask, inc; |
| 1120 | int shift_adj; |
| 1121 | |
| 1122 | a->exp = exp_scalbn(a->exp, scale); |
| 1123 | |
| 1124 | if (a->exp < 0) { |
| 1125 | bool one; |
| 1126 | |
| 1127 | /* All fractional */ |
| 1128 | switch (rmode) { |
| 1129 | case float_round_nearest_even: |
| 1130 | case float_round_nearest_even_max: |
| 1131 | one = false; |
| 1132 | if (a->exp == -1) { |
| 1133 | FloatPartsN tmp; |
| 1134 | /* Shift left one, discarding DECOMPOSED_IMPLICIT_BIT */ |
| 1135 | fracN(add)(&tmp, a, a); |
| 1136 | /* Anything remaining means frac > 0.5. */ |
| 1137 | one = !fracN(eqz)(&tmp); |
| 1138 | } |
| 1139 | break; |
| 1140 | case float_round_ties_away: |
| 1141 | one = a->exp == -1; |
| 1142 | break; |
| 1143 | case float_round_to_zero: |
| 1144 | one = false; |
| 1145 | break; |
| 1146 | case float_round_up: |
| 1147 | one = !a->sign; |
| 1148 | break; |
| 1149 | case float_round_down: |
| 1150 | one = a->sign; |
| 1151 | break; |
| 1152 | case float_round_to_odd: |
| 1153 | case float_round_to_odd_inf: |
| 1154 | one = true; |
| 1155 | break; |
| 1156 | default: |
| 1157 | g_assert_not_reached(); |
| 1158 | } |
| 1159 | |
| 1160 | fracN(clear)(a); |
| 1161 | a->exp = 0; |
| 1162 | if (one) { |
| 1163 | a->frac_hi = DECOMPOSED_IMPLICIT_BIT; |
| 1164 | } else { |
| 1165 | a->cls = float_class_zero; |
| 1166 | } |
| 1167 | return true; |
| 1168 | } |
| 1169 | |
| 1170 | if (N > 64 && a->exp < N - 64) { |
| 1171 | /* |
| 1172 | * Rounding is not in the low word -- shift lsb to bit 2, |
| 1173 | * which leaves room for sticky and rounding bit. |
| 1174 | */ |
| 1175 | shift_adj = (N - 1) - (a->exp + 2); |
| 1176 | fracN(shrjam)(a, shift_adj); |
| 1177 | frac_lsb = 1 << 2; |
| 1178 | } else { |
| 1179 | /* |
| 1180 | * Rounding is in the low word -- compute the lsb offset for rounding |
| 1181 | * and for clamping to the target precision, then map it to an offset |
| 1182 | * within frac_lo. |
| 1183 | */ |
| 1184 | shift_adj = 0; |
| 1185 | frac_lsb = DECOMPOSED_IMPLICIT_BIT >> (MIN(a->exp, frac_size) & 63); |
| 1186 | } |
| 1187 | |
| 1188 | frac_lsbm1 = frac_lsb >> 1; |
| 1189 | rnd_mask = frac_lsb - 1; |
| 1190 | rnd_even_mask = rnd_mask | frac_lsb; |
| 1191 | |
| 1192 | if (!(a->frac_lo & rnd_mask)) { |
| 1193 | /* Fractional bits already clear, undo the shift above. */ |
| 1194 | fracN(shl)(a, shift_adj); |
| 1195 | return false; |
| 1196 | } |
| 1197 | |
| 1198 | switch (rmode) { |
| 1199 | case float_round_nearest_even: |
| 1200 | case float_round_nearest_even_max: |
| 1201 | inc = ((a->frac_lo & rnd_even_mask) != frac_lsbm1 ? frac_lsbm1 : 0); |
| 1202 | break; |
| 1203 | case float_round_ties_away: |
| 1204 | inc = frac_lsbm1; |
| 1205 | break; |
| 1206 | case float_round_to_zero: |
| 1207 | inc = 0; |
| 1208 | break; |
| 1209 | case float_round_up: |
| 1210 | inc = a->sign ? 0 : rnd_mask; |
| 1211 | break; |
| 1212 | case float_round_down: |
| 1213 | inc = a->sign ? rnd_mask : 0; |
| 1214 | break; |
| 1215 | case float_round_to_odd: |
| 1216 | case float_round_to_odd_inf: |
| 1217 | inc = a->frac_lo & frac_lsb ? 0 : rnd_mask; |
| 1218 | break; |
| 1219 | default: |
| 1220 | g_assert_not_reached(); |
| 1221 | } |
| 1222 | |
| 1223 | if (shift_adj == 0) { |
| 1224 | if (fracN(addi)(a, a, inc)) { |
| 1225 | fracN(shr)(a, 1); |
| 1226 | a->frac_hi |= DECOMPOSED_IMPLICIT_BIT; |
| 1227 | a->exp++; |
| 1228 | } |
| 1229 | a->frac_lo &= ~rnd_mask; |
| 1230 | } else { |
| 1231 | fracN(addi)(a, a, inc); |
| 1232 | a->frac_lo &= ~rnd_mask; |
| 1233 | /* Be careful shifting back, not to overflow */ |
| 1234 | fracN(shl)(a, shift_adj - 1); |
| 1235 | if (a->frac_hi & DECOMPOSED_IMPLICIT_BIT) { |
| 1236 | a->exp++; |
| 1237 | } else { |
| 1238 | fracN(add)(a, a, a); |
| 1239 | } |
| 1240 | } |
| 1241 | return true; |
| 1242 | } |
| 1243 | |
| 1244 | FloatPartsN partsN(round_to_int)(const FloatPartsN *a, |
| 1245 | FloatRoundMode rmode, |
| 1246 | int scale, float_status *s, |
| 1247 | const FloatFmt *fmt) |
| 1248 | { |
| 1249 | switch (a->cls) { |
| 1250 | case float_class_qnan: |
| 1251 | case float_class_snan: |
| 1252 | return partsN(return_nan)(a, s); |
| 1253 | case float_class_zero: |
| 1254 | case float_class_inf: |
| 1255 | return *a; |
| 1256 | case float_class_normal: |
| 1257 | case float_class_denormal: |
| 1258 | { |
| 1259 | FloatPartsN r = *a; |
| 1260 | if (partsN(round_to_int_normal)(&r, rmode, scale, fmt->frac_size)) { |
| 1261 | float_raise(float_flag_inexact, s); |
| 1262 | } |
| 1263 | return r; |
| 1264 | } |
| 1265 | default: |
| 1266 | g_assert_not_reached(); |
| 1267 | } |
| 1268 | } |
| 1269 | |
| 1270 | /* |
| 1271 | * Returns the result of converting the floating-point value `a' to |
| 1272 | * the two's complement integer format. The conversion is performed |
| 1273 | * according to the IEC/IEEE Standard for Binary Floating-Point |
| 1274 | * Arithmetic---which means in particular that the conversion is |
| 1275 | * rounded according to the current rounding mode. If `a' is a NaN, |
| 1276 | * the largest positive integer is returned. Otherwise, if the |
| 1277 | * conversion overflows, the largest integer with the same sign as `a' |
| 1278 | * is returned. |
| 1279 | */ |
| 1280 | static int64_t partsN(float_to_sint)(FloatPartsN *p, FloatRoundMode rmode, |
| 1281 | int scale, int64_t min, int64_t max, |
| 1282 | float_status *s) |
| 1283 | { |
| 1284 | FloatExceptionFlags flags = 0; |
| 1285 | uint64_t r; |
| 1286 | |
| 1287 | switch (p->cls) { |
| 1288 | case float_class_snan: |
| 1289 | flags |= float_flag_invalid_snan; |
| 1290 | /* fall through */ |
| 1291 | case float_class_qnan: |
| 1292 | flags |= float_flag_invalid; |
| 1293 | r = max; |
| 1294 | break; |
| 1295 | |
| 1296 | case float_class_inf: |
| 1297 | flags = float_flag_invalid | float_flag_invalid_cvti; |
| 1298 | r = p->sign ? min : max; |
| 1299 | break; |
| 1300 | |
| 1301 | case float_class_zero: |
| 1302 | return 0; |
| 1303 | |
| 1304 | case float_class_normal: |
| 1305 | case float_class_denormal: |
| 1306 | /* TODO: N - 2 is frac_size for rounding; could use input fmt. */ |
| 1307 | if (partsN(round_to_int_normal)(p, rmode, scale, N - 2)) { |
| 1308 | flags = float_flag_inexact; |
| 1309 | } |
| 1310 | |
| 1311 | if (p->exp <= DECOMPOSED_BINARY_POINT) { |
| 1312 | r = p->frac_hi >> (DECOMPOSED_BINARY_POINT - p->exp); |
| 1313 | } else { |
| 1314 | r = UINT64_MAX; |
| 1315 | } |
| 1316 | if (p->sign) { |
| 1317 | if (r <= -(uint64_t)min) { |
| 1318 | r = -r; |
| 1319 | } else { |
| 1320 | flags = float_flag_invalid | float_flag_invalid_cvti; |
| 1321 | r = min; |
| 1322 | } |
| 1323 | } else if (r > max) { |
| 1324 | flags = float_flag_invalid | float_flag_invalid_cvti; |
| 1325 | r = max; |
| 1326 | } |
| 1327 | break; |
| 1328 | |
| 1329 | default: |
| 1330 | g_assert_not_reached(); |
| 1331 | } |
| 1332 | |
| 1333 | float_raise(flags, s); |
| 1334 | return r; |
| 1335 | } |
| 1336 | |
| 1337 | /* |
| 1338 | * Returns the result of converting the floating-point value `a' to |
| 1339 | * the unsigned integer format. The conversion is performed according |
| 1340 | * to the IEC/IEEE Standard for Binary Floating-Point |
| 1341 | * Arithmetic---which means in particular that the conversion is |
| 1342 | * rounded according to the current rounding mode. If `a' is a NaN, |
| 1343 | * the largest unsigned integer is returned. Otherwise, if the |
| 1344 | * conversion overflows, the largest unsigned integer is returned. If |
| 1345 | * the 'a' is negative, the result is rounded and zero is returned; |
| 1346 | * values that do not round to zero will raise the inexact exception |
| 1347 | * flag. |
| 1348 | */ |
| 1349 | static uint64_t partsN(float_to_uint)(FloatPartsN *p, FloatRoundMode rmode, |
| 1350 | int scale, uint64_t max, float_status *s) |
| 1351 | { |
| 1352 | FloatExceptionFlags flags = 0; |
| 1353 | uint64_t r; |
| 1354 | |
| 1355 | switch (p->cls) { |
| 1356 | case float_class_snan: |
| 1357 | flags |= float_flag_invalid_snan; |
| 1358 | /* fall through */ |
| 1359 | case float_class_qnan: |
| 1360 | flags |= float_flag_invalid; |
| 1361 | r = max; |
| 1362 | break; |
| 1363 | |
| 1364 | case float_class_inf: |
| 1365 | flags = float_flag_invalid | float_flag_invalid_cvti; |
| 1366 | r = p->sign ? 0 : max; |
| 1367 | break; |
| 1368 | |
| 1369 | case float_class_zero: |
| 1370 | return 0; |
| 1371 | |
| 1372 | case float_class_normal: |
| 1373 | case float_class_denormal: |
| 1374 | /* TODO: N - 2 is frac_size for rounding; could use input fmt. */ |
| 1375 | if (partsN(round_to_int_normal)(p, rmode, scale, N - 2)) { |
| 1376 | flags = float_flag_inexact; |
| 1377 | if (p->cls == float_class_zero) { |
| 1378 | r = 0; |
| 1379 | break; |
| 1380 | } |
| 1381 | } |
| 1382 | |
| 1383 | if (p->sign) { |
| 1384 | flags = float_flag_invalid | float_flag_invalid_cvti; |
| 1385 | r = 0; |
| 1386 | } else if (p->exp > DECOMPOSED_BINARY_POINT) { |
| 1387 | flags = float_flag_invalid | float_flag_invalid_cvti; |
| 1388 | r = max; |
| 1389 | } else { |
| 1390 | r = p->frac_hi >> (DECOMPOSED_BINARY_POINT - p->exp); |
| 1391 | if (r > max) { |
| 1392 | flags = float_flag_invalid | float_flag_invalid_cvti; |
| 1393 | r = max; |
| 1394 | } |
| 1395 | } |
| 1396 | break; |
| 1397 | |
| 1398 | default: |
| 1399 | g_assert_not_reached(); |
| 1400 | } |
| 1401 | |
| 1402 | float_raise(flags, s); |
| 1403 | return r; |
| 1404 | } |
| 1405 | |
| 1406 | /* |
| 1407 | * Integer to float conversions |
| 1408 | * |
| 1409 | * Returns the result of converting the two's complement integer `a' |
| 1410 | * to the floating-point format. The conversion is performed according |
| 1411 | * to the IEC/IEEE Standard for Binary Floating-Point Arithmetic. |
| 1412 | */ |
| 1413 | static void partsN(sint_to_float)(FloatPartsN *p, int64_t a, |
| 1414 | int scale, float_status *s) |
| 1415 | { |
| 1416 | uint64_t f = a; |
| 1417 | int shift; |
| 1418 | |
| 1419 | memset(p, 0, sizeof(*p)); |
| 1420 | |
| 1421 | if (a == 0) { |
| 1422 | p->cls = float_class_zero; |
| 1423 | return; |
| 1424 | } |
| 1425 | |
| 1426 | p->cls = float_class_normal; |
| 1427 | if (a < 0) { |
| 1428 | f = -f; |
| 1429 | p->sign = true; |
| 1430 | } |
| 1431 | shift = clz64(f); |
| 1432 | scale = MIN(MAX(scale, -0x10000), 0x10000); |
| 1433 | |
| 1434 | p->exp = DECOMPOSED_BINARY_POINT - shift + scale; |
| 1435 | p->frac_hi = f << shift; |
| 1436 | } |
| 1437 | |
| 1438 | /* |
| 1439 | * Unsigned Integer to float conversions |
| 1440 | * |
| 1441 | * Returns the result of converting the unsigned integer `a' to the |
| 1442 | * floating-point format. The conversion is performed according to the |
| 1443 | * IEC/IEEE Standard for Binary Floating-Point Arithmetic. |
| 1444 | */ |
| 1445 | static void partsN(uint_to_float)(FloatPartsN *p, uint64_t a, |
| 1446 | int scale, float_status *status) |
| 1447 | { |
| 1448 | memset(p, 0, sizeof(*p)); |
| 1449 | |
| 1450 | if (a == 0) { |
| 1451 | p->cls = float_class_zero; |
| 1452 | } else { |
| 1453 | int shift = clz64(a); |
| 1454 | scale = MIN(MAX(scale, -0x10000), 0x10000); |
| 1455 | p->cls = float_class_normal; |
| 1456 | p->exp = DECOMPOSED_BINARY_POINT - shift + scale; |
| 1457 | p->frac_hi = a << shift; |
| 1458 | } |
| 1459 | } |
| 1460 | |
| 1461 | /* |
| 1462 | * Float min/max. |
| 1463 | */ |
| 1464 | static FloatPartsN *partsN(minmax)(FloatPartsN *a, FloatPartsN *b, |
| 1465 | float_status *s, int flags) |
| 1466 | { |
| 1467 | int ab_mask = float_cmask(a->cls) | float_cmask(b->cls); |
| 1468 | int a_exp, b_exp, cmp; |
| 1469 | |
| 1470 | if (unlikely(ab_mask & float_cmask_anynan)) { |
| 1471 | /* |
| 1472 | * For minNum/maxNum (IEEE 754-2008) |
| 1473 | * or minimumNumber/maximumNumber (IEEE 754-2019), |
| 1474 | * if one operand is a QNaN, and the other |
| 1475 | * operand is numerical, then return numerical argument. |
| 1476 | */ |
| 1477 | if ((flags & (float_minmax_isnum | float_minmax_isnumber)) |
| 1478 | && !(ab_mask & float_cmask_snan) |
| 1479 | && (ab_mask & ~float_cmask_qnan)) { |
| 1480 | record_denormals_used(ab_mask, s); |
| 1481 | return is_nan(a->cls) ? b : a; |
| 1482 | } |
| 1483 | |
| 1484 | /* |
| 1485 | * In IEEE 754-2019, minNum, maxNum, minNumMag and maxNumMag |
| 1486 | * are removed and replaced with minimum, minimumNumber, maximum |
| 1487 | * and maximumNumber. |
| 1488 | * minimumNumber/maximumNumber behavior for SNaN is changed to: |
| 1489 | * If both operands are NaNs, a QNaN is returned. |
| 1490 | * If either operand is a SNaN, |
| 1491 | * an invalid operation exception is signaled, |
| 1492 | * but unless both operands are NaNs, |
| 1493 | * the SNaN is otherwise ignored and not converted to a QNaN. |
| 1494 | */ |
| 1495 | if ((flags & float_minmax_isnumber) |
| 1496 | && (ab_mask & float_cmask_snan) |
| 1497 | && (ab_mask & ~float_cmask_anynan)) { |
| 1498 | float_raise(float_flag_invalid, s); |
| 1499 | return is_nan(a->cls) ? b : a; |
| 1500 | } |
| 1501 | |
| 1502 | *a = partsN(pick_nan)(a, b, s); |
| 1503 | return a; |
| 1504 | } |
| 1505 | |
| 1506 | record_denormals_used(ab_mask, s); |
| 1507 | |
| 1508 | a_exp = a->exp; |
| 1509 | b_exp = b->exp; |
| 1510 | |
| 1511 | if (unlikely(!cmask_is_only_normals(ab_mask))) { |
| 1512 | switch (a->cls) { |
| 1513 | case float_class_normal: |
| 1514 | case float_class_denormal: |
| 1515 | break; |
| 1516 | case float_class_inf: |
| 1517 | a_exp = INT16_MAX; |
| 1518 | break; |
| 1519 | case float_class_zero: |
| 1520 | a_exp = INT16_MIN; |
| 1521 | break; |
| 1522 | default: |
| 1523 | g_assert_not_reached(); |
| 1524 | } |
| 1525 | switch (b->cls) { |
| 1526 | case float_class_normal: |
| 1527 | case float_class_denormal: |
| 1528 | break; |
| 1529 | case float_class_inf: |
| 1530 | b_exp = INT16_MAX; |
| 1531 | break; |
| 1532 | case float_class_zero: |
| 1533 | b_exp = INT16_MIN; |
| 1534 | break; |
| 1535 | default: |
| 1536 | g_assert_not_reached(); |
| 1537 | } |
| 1538 | } |
| 1539 | |
| 1540 | /* Compare magnitudes. */ |
| 1541 | cmp = a_exp - b_exp; |
| 1542 | if (cmp == 0) { |
| 1543 | cmp = fracN(cmp)(a, b); |
| 1544 | } |
| 1545 | |
| 1546 | /* |
| 1547 | * Take the sign into account. |
| 1548 | * For ismag, only do this if the magnitudes are equal. |
| 1549 | */ |
| 1550 | if (!(flags & float_minmax_ismag) || cmp == 0) { |
| 1551 | if (a->sign != b->sign) { |
| 1552 | /* For differing signs, the negative operand is less. */ |
| 1553 | cmp = a->sign ? -1 : 1; |
| 1554 | } else if (a->sign) { |
| 1555 | /* For two negative operands, invert the magnitude comparison. */ |
| 1556 | cmp = -cmp; |
| 1557 | } |
| 1558 | } |
| 1559 | |
| 1560 | if (flags & float_minmax_ismin) { |
| 1561 | cmp = -cmp; |
| 1562 | } |
| 1563 | return cmp < 0 ? b : a; |
| 1564 | } |
| 1565 | |
| 1566 | /* |
| 1567 | * Floating point compare |
| 1568 | */ |
| 1569 | FloatRelation partsN(compare)(const FloatPartsN *a, const FloatPartsN *b, |
| 1570 | float_status *s, bool is_quiet) |
| 1571 | { |
| 1572 | int ab_mask = float_cmask(a->cls) | float_cmask(b->cls); |
| 1573 | |
| 1574 | if (likely(cmask_is_only_normals(ab_mask))) { |
| 1575 | FloatRelation cmp; |
| 1576 | |
| 1577 | record_denormals_used(ab_mask, s); |
| 1578 | |
| 1579 | if (a->sign != b->sign) { |
| 1580 | goto a_sign; |
| 1581 | } |
| 1582 | if (a->exp == b->exp) { |
| 1583 | cmp = fracN(cmp)(a, b); |
| 1584 | } else if (a->exp < b->exp) { |
| 1585 | cmp = float_relation_less; |
| 1586 | } else { |
| 1587 | cmp = float_relation_greater; |
| 1588 | } |
| 1589 | if (a->sign) { |
| 1590 | cmp = -cmp; |
| 1591 | } |
| 1592 | return cmp; |
| 1593 | } |
| 1594 | |
| 1595 | if (unlikely(ab_mask & float_cmask_anynan)) { |
| 1596 | if (ab_mask & float_cmask_snan) { |
| 1597 | float_raise(float_flag_invalid | float_flag_invalid_snan, s); |
| 1598 | } else if (!is_quiet) { |
| 1599 | float_raise(float_flag_invalid, s); |
| 1600 | } |
| 1601 | return float_relation_unordered; |
| 1602 | } |
| 1603 | |
| 1604 | record_denormals_used(ab_mask, s); |
| 1605 | |
| 1606 | if (ab_mask & float_cmask_zero) { |
| 1607 | if (ab_mask == float_cmask_zero) { |
| 1608 | return float_relation_equal; |
| 1609 | } else if (a->cls == float_class_zero) { |
| 1610 | goto b_sign; |
| 1611 | } else { |
| 1612 | goto a_sign; |
| 1613 | } |
| 1614 | } |
| 1615 | |
| 1616 | if (ab_mask == float_cmask_inf) { |
| 1617 | if (a->sign == b->sign) { |
| 1618 | return float_relation_equal; |
| 1619 | } |
| 1620 | } else if (b->cls == float_class_inf) { |
| 1621 | goto b_sign; |
| 1622 | } else { |
| 1623 | g_assert(a->cls == float_class_inf); |
| 1624 | } |
| 1625 | |
| 1626 | a_sign: |
| 1627 | return a->sign ? float_relation_less : float_relation_greater; |
| 1628 | b_sign: |
| 1629 | return b->sign ? float_relation_greater : float_relation_less; |
| 1630 | } |
| 1631 | |
| 1632 | /* |
| 1633 | * Multiply A by 2 raised to the power N. |
| 1634 | */ |
| 1635 | FloatPartsN partsN(scalbn)(const FloatPartsN *a, int n, float_status *s) |
| 1636 | { |
| 1637 | switch (a->cls) { |
| 1638 | case float_class_snan: |
| 1639 | case float_class_qnan: |
| 1640 | return partsN(return_nan)(a, s); |
| 1641 | case float_class_zero: |
| 1642 | case float_class_inf: |
| 1643 | return *a; |
| 1644 | case float_class_denormal: |
| 1645 | float_raise(float_flag_input_denormal_used, s); |
| 1646 | /* fall through */ |
| 1647 | case float_class_normal: |
| 1648 | { |
| 1649 | FloatPartsN r = *a; |
| 1650 | r.exp = exp_scalbn(r.exp, n); |
| 1651 | return r; |
| 1652 | } |
| 1653 | default: |
| 1654 | g_assert_not_reached(); |
| 1655 | } |
| 1656 | } |