qemu

softfloat-parts.c.inc (44581B)

---

 /*
  * QEMU float support
  *
  * The code in this source file is derived from release 2a of the SoftFloat
  * IEC/IEEE Floating-point Arithmetic Package. Those parts of the code (and
  * some later contributions) are provided under that license, as detailed below.
  * It has subsequently been modified by contributors to the QEMU Project,
  * so some portions are provided under:
  *  the SoftFloat-2a license
  *  the BSD license
  *  GPL-v2-or-later
  *
  * Any future contributions to this file after December 1st 2014 will be
  * taken to be licensed under the Softfloat-2a license unless specifically
  * indicated otherwise.
  */
 
 static void partsN(return_nan)(FloatPartsN *a, float_status *s)
 {
     switch (a->cls) {
     case float_class_snan:
         float_raise(float_flag_invalid | float_flag_invalid_snan, s);
         if (s->default_nan_mode) {
             parts_default_nan(a, s);
         } else {
             parts_silence_nan(a, s);
         }
         break;
     case float_class_qnan:
         if (s->default_nan_mode) {
             parts_default_nan(a, s);
         }
         break;
     default:
         g_assert_not_reached();
     }
 }
 
 static FloatPartsN *partsN(pick_nan)(FloatPartsN *a, FloatPartsN *b,
                                      float_status *s)
 {
     if (is_snan(a->cls) || is_snan(b->cls)) {
         float_raise(float_flag_invalid | float_flag_invalid_snan, s);
     }
 
     if (s->default_nan_mode) {
         parts_default_nan(a, s);
     } else {
         int cmp = frac_cmp(a, b);
         if (cmp == 0) {
             cmp = a->sign < b->sign;
         }
 
         if (pickNaN(a->cls, b->cls, cmp > 0, s)) {
             a = b;
         }
         if (is_snan(a->cls)) {
             parts_silence_nan(a, s);
         }
     }
     return a;
 }
 
 static FloatPartsN *partsN(pick_nan_muladd)(FloatPartsN *a, FloatPartsN *b,
                                             FloatPartsN *c, float_status *s,
                                             int ab_mask, int abc_mask)
 {
     int which;
 
     if (unlikely(abc_mask & float_cmask_snan)) {
         float_raise(float_flag_invalid | float_flag_invalid_snan, s);
     }
 
     which = pickNaNMulAdd(a->cls, b->cls, c->cls,
                           ab_mask == float_cmask_infzero, s);
 
     if (s->default_nan_mode || which == 3) {
         /*
          * Note that this check is after pickNaNMulAdd so that function
          * has an opportunity to set the Invalid flag for infzero.
          */
         parts_default_nan(a, s);
         return a;
     }
 
     switch (which) {
     case 0:
         break;
     case 1:
         a = b;
         break;
     case 2:
         a = c;
         break;
     default:
         g_assert_not_reached();
     }
     if (is_snan(a->cls)) {
         parts_silence_nan(a, s);
     }
     return a;
 }
 
 /*
  * Canonicalize the FloatParts structure.  Determine the class,
  * unbias the exponent, and normalize the fraction.
  */
 static void partsN(canonicalize)(FloatPartsN *p, float_status *status,
                                  const FloatFmt *fmt)
 {
     if (unlikely(p->exp == 0)) {
         if (likely(frac_eqz(p))) {
             p->cls = float_class_zero;
         } else if (status->flush_inputs_to_zero) {
             float_raise(float_flag_input_denormal, status);
             p->cls = float_class_zero;
             frac_clear(p);
         } else {
             int shift = frac_normalize(p);
             p->cls = float_class_normal;
             p->exp = fmt->frac_shift - fmt->exp_bias - shift + 1;
         }
     } else if (likely(p->exp < fmt->exp_max) || fmt->arm_althp) {
         p->cls = float_class_normal;
         p->exp -= fmt->exp_bias;
         frac_shl(p, fmt->frac_shift);
         p->frac_hi |= DECOMPOSED_IMPLICIT_BIT;
     } else if (likely(frac_eqz(p))) {
         p->cls = float_class_inf;
     } else {
         frac_shl(p, fmt->frac_shift);
         p->cls = (parts_is_snan_frac(p->frac_hi, status)
                   ? float_class_snan : float_class_qnan);
     }
 }
 
 /*
  * Round and uncanonicalize a floating-point number by parts. There
  * are FRAC_SHIFT bits that may require rounding at the bottom of the
  * fraction; these bits will be removed. The exponent will be biased
  * by EXP_BIAS and must be bounded by [EXP_MAX-1, 0].
  */
 static void partsN(uncanon_normal)(FloatPartsN *p, float_status *s,
                                    const FloatFmt *fmt)
 {
     const int exp_max = fmt->exp_max;
     const int frac_shift = fmt->frac_shift;
     const uint64_t round_mask = fmt->round_mask;
     const uint64_t frac_lsb = round_mask + 1;
     const uint64_t frac_lsbm1 = round_mask ^ (round_mask >> 1);
     const uint64_t roundeven_mask = round_mask | frac_lsb;
     uint64_t inc;
     bool overflow_norm = false;
     int exp, flags = 0;
 
     switch (s->float_rounding_mode) {
     case float_round_nearest_even:
         if (N > 64 && frac_lsb == 0) {
             inc = ((p->frac_hi & 1) || (p->frac_lo & round_mask) != frac_lsbm1
                    ? frac_lsbm1 : 0);
         } else {
             inc = ((p->frac_lo & roundeven_mask) != frac_lsbm1
                    ? frac_lsbm1 : 0);
         }
         break;
     case float_round_ties_away:
         inc = frac_lsbm1;
         break;
     case float_round_to_zero:
         overflow_norm = true;
         inc = 0;
         break;
     case float_round_up:
         inc = p->sign ? 0 : round_mask;
         overflow_norm = p->sign;
         break;
     case float_round_down:
         inc = p->sign ? round_mask : 0;
         overflow_norm = !p->sign;
         break;
     case float_round_to_odd:
         overflow_norm = true;
         /* fall through */
     case float_round_to_odd_inf:
         if (N > 64 && frac_lsb == 0) {
             inc = p->frac_hi & 1 ? 0 : round_mask;
         } else {
             inc = p->frac_lo & frac_lsb ? 0 : round_mask;
         }
         break;
     default:
         g_assert_not_reached();
     }
 
     exp = p->exp + fmt->exp_bias;
     if (likely(exp > 0)) {
         if (p->frac_lo & round_mask) {
             flags |= float_flag_inexact;
             if (frac_addi(p, p, inc)) {
                 frac_shr(p, 1);
                 p->frac_hi |= DECOMPOSED_IMPLICIT_BIT;
                 exp++;
             }
             p->frac_lo &= ~round_mask;
         }
 
         if (fmt->arm_althp) {
             /* ARM Alt HP eschews Inf and NaN for a wider exponent.  */
             if (unlikely(exp > exp_max)) {
                 /* Overflow.  Return the maximum normal.  */
                 flags = float_flag_invalid;
                 exp = exp_max;
                 frac_allones(p);
                 p->frac_lo &= ~round_mask;
             }
         } else if (unlikely(exp >= exp_max)) {
             flags |= float_flag_overflow;
             if (s->rebias_overflow) {
                 exp -= fmt->exp_re_bias;
             } else if (overflow_norm) {
                 flags |= float_flag_inexact;
                 exp = exp_max - 1;
                 frac_allones(p);
                 p->frac_lo &= ~round_mask;
             } else {
                 flags |= float_flag_inexact;
                 p->cls = float_class_inf;
                 exp = exp_max;
                 frac_clear(p);
             }
         }
         frac_shr(p, frac_shift);
     } else if (unlikely(s->rebias_underflow)) {
         flags |= float_flag_underflow;
         exp += fmt->exp_re_bias;
         if (p->frac_lo & round_mask) {
             flags |= float_flag_inexact;
             if (frac_addi(p, p, inc)) {
                 frac_shr(p, 1);
                 p->frac_hi |= DECOMPOSED_IMPLICIT_BIT;
                 exp++;
             }
             p->frac_lo &= ~round_mask;
         }
         frac_shr(p, frac_shift);
     } else if (s->flush_to_zero) {
         flags |= float_flag_output_denormal;
         p->cls = float_class_zero;
         exp = 0;
         frac_clear(p);
     } else {
         bool is_tiny = s->tininess_before_rounding || exp < 0;
 
         if (!is_tiny) {
             FloatPartsN discard;
             is_tiny = !frac_addi(&discard, p, inc);
         }
 
         frac_shrjam(p, 1 - exp);
 
         if (p->frac_lo & round_mask) {
             /* Need to recompute round-to-even/round-to-odd. */
             switch (s->float_rounding_mode) {
             case float_round_nearest_even:
                 if (N > 64 && frac_lsb == 0) {
                     inc = ((p->frac_hi & 1) ||
                            (p->frac_lo & round_mask) != frac_lsbm1
                            ? frac_lsbm1 : 0);
                 } else {
                     inc = ((p->frac_lo & roundeven_mask) != frac_lsbm1
                            ? frac_lsbm1 : 0);
                 }
                 break;
             case float_round_to_odd:
             case float_round_to_odd_inf:
                 if (N > 64 && frac_lsb == 0) {
                     inc = p->frac_hi & 1 ? 0 : round_mask;
                 } else {
                     inc = p->frac_lo & frac_lsb ? 0 : round_mask;
                 }
                 break;
             default:
                 break;
             }
             flags |= float_flag_inexact;
             frac_addi(p, p, inc);
             p->frac_lo &= ~round_mask;
         }
 
         exp = (p->frac_hi & DECOMPOSED_IMPLICIT_BIT) != 0;
         frac_shr(p, frac_shift);
 
         if (is_tiny && (flags & float_flag_inexact)) {
             flags |= float_flag_underflow;
         }
         if (exp == 0 && frac_eqz(p)) {
             p->cls = float_class_zero;
         }
     }
     p->exp = exp;
     float_raise(flags, s);
 }
 
 static void partsN(uncanon)(FloatPartsN *p, float_status *s,
                             const FloatFmt *fmt)
 {
     if (likely(p->cls == float_class_normal)) {
         parts_uncanon_normal(p, s, fmt);
     } else {
         switch (p->cls) {
         case float_class_zero:
             p->exp = 0;
             frac_clear(p);
             return;
         case float_class_inf:
             g_assert(!fmt->arm_althp);
             p->exp = fmt->exp_max;
             frac_clear(p);
             return;
         case float_class_qnan:
         case float_class_snan:
             g_assert(!fmt->arm_althp);
             p->exp = fmt->exp_max;
             frac_shr(p, fmt->frac_shift);
             return;
         default:
             break;
         }
         g_assert_not_reached();
     }
 }
 
 /*
  * Returns the result of adding or subtracting the values of the
  * floating-point values `a' and `b'. The operation is performed
  * according to the IEC/IEEE Standard for Binary Floating-Point
  * Arithmetic.
  */
 static FloatPartsN *partsN(addsub)(FloatPartsN *a, FloatPartsN *b,
                                    float_status *s, bool subtract)
 {
     bool b_sign = b->sign ^ subtract;
     int ab_mask = float_cmask(a->cls) | float_cmask(b->cls);
 
     if (a->sign != b_sign) {
         /* Subtraction */
         if (likely(ab_mask == float_cmask_normal)) {
             if (parts_sub_normal(a, b)) {
                 return a;
             }
             /* Subtract was exact, fall through to set sign. */
             ab_mask = float_cmask_zero;
         }
 
         if (ab_mask == float_cmask_zero) {
             a->sign = s->float_rounding_mode == float_round_down;
             return a;
         }
 
         if (unlikely(ab_mask & float_cmask_anynan)) {
             goto p_nan;
         }
 
         if (ab_mask & float_cmask_inf) {
             if (a->cls != float_class_inf) {
                 /* N - Inf */
                 goto return_b;
             }
             if (b->cls != float_class_inf) {
                 /* Inf - N */
                 return a;
             }
             /* Inf - Inf */
             float_raise(float_flag_invalid | float_flag_invalid_isi, s);
             parts_default_nan(a, s);
             return a;
         }
     } else {
         /* Addition */
         if (likely(ab_mask == float_cmask_normal)) {
             parts_add_normal(a, b);
             return a;
         }
 
         if (ab_mask == float_cmask_zero) {
             return a;
         }
 
         if (unlikely(ab_mask & float_cmask_anynan)) {
             goto p_nan;
         }
 
         if (ab_mask & float_cmask_inf) {
             a->cls = float_class_inf;
             return a;
         }
     }
 
     if (b->cls == float_class_zero) {
         g_assert(a->cls == float_class_normal);
         return a;
     }
 
     g_assert(a->cls == float_class_zero);
     g_assert(b->cls == float_class_normal);
  return_b:
     b->sign = b_sign;
     return b;
 
  p_nan:
     return parts_pick_nan(a, b, s);
 }
 
 /*
  * Returns the result of multiplying the floating-point values `a' and
  * `b'. The operation is performed according to the IEC/IEEE Standard
  * for Binary Floating-Point Arithmetic.
  */
 static FloatPartsN *partsN(mul)(FloatPartsN *a, FloatPartsN *b,
                                 float_status *s)
 {
     int ab_mask = float_cmask(a->cls) | float_cmask(b->cls);
     bool sign = a->sign ^ b->sign;
 
     if (likely(ab_mask == float_cmask_normal)) {
         FloatPartsW tmp;
 
         frac_mulw(&tmp, a, b);
         frac_truncjam(a, &tmp);
 
         a->exp += b->exp + 1;
         if (!(a->frac_hi & DECOMPOSED_IMPLICIT_BIT)) {
             frac_add(a, a, a);
             a->exp -= 1;
         }
 
         a->sign = sign;
         return a;
     }
 
     /* Inf * Zero == NaN */
     if (unlikely(ab_mask == float_cmask_infzero)) {
         float_raise(float_flag_invalid | float_flag_invalid_imz, s);
         parts_default_nan(a, s);
         return a;
     }
 
     if (unlikely(ab_mask & float_cmask_anynan)) {
         return parts_pick_nan(a, b, s);
     }
 
     /* Multiply by 0 or Inf */
     if (ab_mask & float_cmask_inf) {
         a->cls = float_class_inf;
         a->sign = sign;
         return a;
     }
 
     g_assert(ab_mask & float_cmask_zero);
     a->cls = float_class_zero;
     a->sign = sign;
     return a;
 }
 
 /*
  * Returns the result of multiplying the floating-point values `a' and
  * `b' then adding 'c', with no intermediate rounding step after the
  * multiplication. The operation is performed according to the
  * IEC/IEEE Standard for Binary Floating-Point Arithmetic 754-2008.
  * The flags argument allows the caller to select negation of the
  * addend, the intermediate product, or the final result. (The
  * difference between this and having the caller do a separate
  * negation is that negating externally will flip the sign bit on NaNs.)
  *
  * Requires A and C extracted into a double-sized structure to provide the
  * extra space for the widening multiply.
  */
 static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b,
                                    FloatPartsN *c, int flags, float_status *s)
 {
     int ab_mask, abc_mask;
     FloatPartsW p_widen, c_widen;
 
     ab_mask = float_cmask(a->cls) | float_cmask(b->cls);
     abc_mask = float_cmask(c->cls) | ab_mask;
 
     /*
      * It is implementation-defined whether the cases of (0,inf,qnan)
      * and (inf,0,qnan) raise InvalidOperation or not (and what QNaN
      * they return if they do), so we have to hand this information
      * off to the target-specific pick-a-NaN routine.
      */
     if (unlikely(abc_mask & float_cmask_anynan)) {
         return parts_pick_nan_muladd(a, b, c, s, ab_mask, abc_mask);
     }
 
     if (flags & float_muladd_negate_c) {
         c->sign ^= 1;
     }
 
     /* Compute the sign of the product into A. */
     a->sign ^= b->sign;
     if (flags & float_muladd_negate_product) {
         a->sign ^= 1;
     }
 
     if (unlikely(ab_mask != float_cmask_normal)) {
         if (unlikely(ab_mask == float_cmask_infzero)) {
             float_raise(float_flag_invalid | float_flag_invalid_imz, s);
             goto d_nan;
         }
 
         if (ab_mask & float_cmask_inf) {
             if (c->cls == float_class_inf && a->sign != c->sign) {
                 float_raise(float_flag_invalid | float_flag_invalid_isi, s);
                 goto d_nan;
             }
             goto return_inf;
         }
 
         g_assert(ab_mask & float_cmask_zero);
         if (c->cls == float_class_normal) {
             *a = *c;
             goto return_normal;
         }
         if (c->cls == float_class_zero) {
             if (a->sign != c->sign) {
                 goto return_sub_zero;
             }
             goto return_zero;
         }
         g_assert(c->cls == float_class_inf);
     }
 
     if (unlikely(c->cls == float_class_inf)) {
         a->sign = c->sign;
         goto return_inf;
     }
 
     /* Perform the multiplication step. */
     p_widen.sign = a->sign;
     p_widen.exp = a->exp + b->exp + 1;
     frac_mulw(&p_widen, a, b);
     if (!(p_widen.frac_hi & DECOMPOSED_IMPLICIT_BIT)) {
         frac_add(&p_widen, &p_widen, &p_widen);
         p_widen.exp -= 1;
     }
 
     /* Perform the addition step. */
     if (c->cls != float_class_zero) {
         /* Zero-extend C to less significant bits. */
         frac_widen(&c_widen, c);
         c_widen.exp = c->exp;
 
         if (a->sign == c->sign) {
             parts_add_normal(&p_widen, &c_widen);
         } else if (!parts_sub_normal(&p_widen, &c_widen)) {
             goto return_sub_zero;
         }
     }
 
     /* Narrow with sticky bit, for proper rounding later. */
     frac_truncjam(a, &p_widen);
     a->sign = p_widen.sign;
     a->exp = p_widen.exp;
 
  return_normal:
     if (flags & float_muladd_halve_result) {
         a->exp -= 1;
     }
  finish_sign:
     if (flags & float_muladd_negate_result) {
         a->sign ^= 1;
     }
     return a;
 
  return_sub_zero:
     a->sign = s->float_rounding_mode == float_round_down;
  return_zero:
     a->cls = float_class_zero;
     goto finish_sign;
 
  return_inf:
     a->cls = float_class_inf;
     goto finish_sign;
 
  d_nan:
     parts_default_nan(a, s);
     return a;
 }
 
 /*
  * Returns the result of dividing the floating-point value `a' by the
  * corresponding value `b'. The operation is performed according to
  * the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
  */
 static FloatPartsN *partsN(div)(FloatPartsN *a, FloatPartsN *b,
                                 float_status *s)
 {
     int ab_mask = float_cmask(a->cls) | float_cmask(b->cls);
     bool sign = a->sign ^ b->sign;
 
     if (likely(ab_mask == float_cmask_normal)) {
         a->sign = sign;
         a->exp -= b->exp + frac_div(a, b);
         return a;
     }
 
     /* 0/0 or Inf/Inf => NaN */
     if (unlikely(ab_mask == float_cmask_zero)) {
         float_raise(float_flag_invalid | float_flag_invalid_zdz, s);
         goto d_nan;
     }
     if (unlikely(ab_mask == float_cmask_inf)) {
         float_raise(float_flag_invalid | float_flag_invalid_idi, s);
         goto d_nan;
     }
 
     /* All the NaN cases */
     if (unlikely(ab_mask & float_cmask_anynan)) {
         return parts_pick_nan(a, b, s);
     }
 
     a->sign = sign;
 
     /* Inf / X */
     if (a->cls == float_class_inf) {
         return a;
     }
 
     /* 0 / X */
     if (a->cls == float_class_zero) {
         return a;
     }
 
     /* X / Inf */
     if (b->cls == float_class_inf) {
         a->cls = float_class_zero;
         return a;
     }
 
     /* X / 0 => Inf */
     g_assert(b->cls == float_class_zero);
     float_raise(float_flag_divbyzero, s);
     a->cls = float_class_inf;
     return a;
 
  d_nan:
     parts_default_nan(a, s);
     return a;
 }
 
 /*
  * Floating point remainder, per IEC/IEEE, or modulus.
  */
 static FloatPartsN *partsN(modrem)(FloatPartsN *a, FloatPartsN *b,
                                    uint64_t *mod_quot, float_status *s)
 {
     int ab_mask = float_cmask(a->cls) | float_cmask(b->cls);
 
     if (likely(ab_mask == float_cmask_normal)) {
         frac_modrem(a, b, mod_quot);
         return a;
     }
 
     if (mod_quot) {
         *mod_quot = 0;
     }
 
     /* All the NaN cases */
     if (unlikely(ab_mask & float_cmask_anynan)) {
         return parts_pick_nan(a, b, s);
     }
 
     /* Inf % N; N % 0 */
     if (a->cls == float_class_inf || b->cls == float_class_zero) {
         float_raise(float_flag_invalid, s);
         parts_default_nan(a, s);
         return a;
     }
 
     /* N % Inf; 0 % N */
     g_assert(b->cls == float_class_inf || a->cls == float_class_zero);
     return a;
 }
 
 /*
  * Square Root
  *
  * The base algorithm is lifted from
  * https://git.musl-libc.org/cgit/musl/tree/src/math/sqrtf.c
  * https://git.musl-libc.org/cgit/musl/tree/src/math/sqrt.c
  * https://git.musl-libc.org/cgit/musl/tree/src/math/sqrtl.c
  * and is thus MIT licenced.
  */
 static void partsN(sqrt)(FloatPartsN *a, float_status *status,
                          const FloatFmt *fmt)
 {
     const uint32_t three32 = 3u << 30;
     const uint64_t three64 = 3ull << 62;
     uint32_t d32, m32, r32, s32, u32;            /* 32-bit computation */
     uint64_t d64, m64, r64, s64, u64;            /* 64-bit computation */
     uint64_t dh, dl, rh, rl, sh, sl, uh, ul;     /* 128-bit computation */
     uint64_t d0h, d0l, d1h, d1l, d2h, d2l;
     uint64_t discard;
     bool exp_odd;
     size_t index;
 
     if (unlikely(a->cls != float_class_normal)) {
         switch (a->cls) {
         case float_class_snan:
         case float_class_qnan:
             parts_return_nan(a, status);
             return;
         case float_class_zero:
             return;
         case float_class_inf:
             if (unlikely(a->sign)) {
                 goto d_nan;
             }
             return;
         default:
             g_assert_not_reached();
         }
     }
 
     if (unlikely(a->sign)) {
         goto d_nan;
     }
 
     /*
      * Argument reduction.
      * x = 4^e frac; with integer e, and frac in [1, 4)
      * m = frac fixed point at bit 62, since we're in base 4.
      * If base-2 exponent is odd, exchange that for multiply by 2,
      * which results in no shift.
      */
     exp_odd = a->exp & 1;
     index = extract64(a->frac_hi, 57, 6) | (!exp_odd << 6);
     if (!exp_odd) {
         frac_shr(a, 1);
     }
 
     /*
      * Approximate r ~= 1/sqrt(m) and s ~= sqrt(m) when m in [1, 4).
      *
      * Initial estimate:
      * 7-bit lookup table (1-bit exponent and 6-bit significand).
      *
      * The relative error (e = r0*sqrt(m)-1) of a linear estimate
      * (r0 = a*m + b) is |e| < 0.085955 ~ 0x1.6p-4 at best;
      * a table lookup is faster and needs one less iteration.
      * The 7-bit table gives |e| < 0x1.fdp-9.
      *
      * A Newton-Raphson iteration for r is
      *   s = m*r
      *   d = s*r
      *   u = 3 - d
      *   r = r*u/2
      *
      * Fixed point representations:
      *   m, s, d, u, three are all 2.30; r is 0.32
      */
     m64 = a->frac_hi;
     m32 = m64 >> 32;
 
     r32 = rsqrt_tab[index] << 16;
     /* |r*sqrt(m) - 1| < 0x1.FDp-9 */
 
     s32 = ((uint64_t)m32 * r32) >> 32;
     d32 = ((uint64_t)s32 * r32) >> 32;
     u32 = three32 - d32;
 
     if (N == 64) {
         /* float64 or smaller */
 
         r32 = ((uint64_t)r32 * u32) >> 31;
         /* |r*sqrt(m) - 1| < 0x1.7Bp-16 */
 
         s32 = ((uint64_t)m32 * r32) >> 32;
         d32 = ((uint64_t)s32 * r32) >> 32;
         u32 = three32 - d32;
 
         if (fmt->frac_size <= 23) {
             /* float32 or smaller */
 
             s32 = ((uint64_t)s32 * u32) >> 32;  /* 3.29 */
             s32 = (s32 - 1) >> 6;               /* 9.23 */
             /* s < sqrt(m) < s + 0x1.08p-23 */
 
             /* compute nearest rounded result to 2.23 bits */
             uint32_t d0 = (m32 << 16) - s32 * s32;
             uint32_t d1 = s32 - d0;
             uint32_t d2 = d1 + s32 + 1;
             s32 += d1 >> 31;
             a->frac_hi = (uint64_t)s32 << (64 - 25);
 
             /* increment or decrement for inexact */
             if (d2 != 0) {
                 a->frac_hi += ((int32_t)(d1 ^ d2) < 0 ? -1 : 1);
             }
             goto done;
         }
 
         /* float64 */
 
         r64 = (uint64_t)r32 * u32 * 2;
         /* |r*sqrt(m) - 1| < 0x1.37-p29; convert to 64-bit arithmetic */
         mul64To128(m64, r64, &s64, &discard);
         mul64To128(s64, r64, &d64, &discard);
         u64 = three64 - d64;
 
         mul64To128(s64, u64, &s64, &discard);  /* 3.61 */
         s64 = (s64 - 2) >> 9;                  /* 12.52 */
 
         /* Compute nearest rounded result */
         uint64_t d0 = (m64 << 42) - s64 * s64;
         uint64_t d1 = s64 - d0;
         uint64_t d2 = d1 + s64 + 1;
         s64 += d1 >> 63;
         a->frac_hi = s64 << (64 - 54);
 
         /* increment or decrement for inexact */
         if (d2 != 0) {
             a->frac_hi += ((int64_t)(d1 ^ d2) < 0 ? -1 : 1);
         }
         goto done;
     }
 
     r64 = (uint64_t)r32 * u32 * 2;
     /* |r*sqrt(m) - 1| < 0x1.7Bp-16; convert to 64-bit arithmetic */
 
     mul64To128(m64, r64, &s64, &discard);
     mul64To128(s64, r64, &d64, &discard);
     u64 = three64 - d64;
     mul64To128(u64, r64, &r64, &discard);
     r64 <<= 1;
     /* |r*sqrt(m) - 1| < 0x1.a5p-31 */
 
     mul64To128(m64, r64, &s64, &discard);
     mul64To128(s64, r64, &d64, &discard);
     u64 = three64 - d64;
     mul64To128(u64, r64, &rh, &rl);
     add128(rh, rl, rh, rl, &rh, &rl);
     /* |r*sqrt(m) - 1| < 0x1.c001p-59; change to 128-bit arithmetic */
 
     mul128To256(a->frac_hi, a->frac_lo, rh, rl, &sh, &sl, &discard, &discard);
     mul128To256(sh, sl, rh, rl, &dh, &dl, &discard, &discard);
     sub128(three64, 0, dh, dl, &uh, &ul);
     mul128To256(uh, ul, sh, sl, &sh, &sl, &discard, &discard);  /* 3.125 */
     /* -0x1p-116 < s - sqrt(m) < 0x3.8001p-125 */
 
     sub128(sh, sl, 0, 4, &sh, &sl);
     shift128Right(sh, sl, 13, &sh, &sl);  /* 16.112 */
     /* s < sqrt(m) < s + 1ulp */
 
     /* Compute nearest rounded result */
     mul64To128(sl, sl, &d0h, &d0l);
     d0h += 2 * sh * sl;
     sub128(a->frac_lo << 34, 0, d0h, d0l, &d0h, &d0l);
     sub128(sh, sl, d0h, d0l, &d1h, &d1l);
     add128(sh, sl, 0, 1, &d2h, &d2l);
     add128(d2h, d2l, d1h, d1l, &d2h, &d2l);
     add128(sh, sl, 0, d1h >> 63, &sh, &sl);
     shift128Left(sh, sl, 128 - 114, &sh, &sl);
 
     /* increment or decrement for inexact */
     if (d2h | d2l) {
         if ((int64_t)(d1h ^ d2h) < 0) {
             sub128(sh, sl, 0, 1, &sh, &sl);
         } else {
             add128(sh, sl, 0, 1, &sh, &sl);
         }
     }
     a->frac_lo = sl;
     a->frac_hi = sh;
 
  done:
     /* Convert back from base 4 to base 2. */
     a->exp >>= 1;
     if (!(a->frac_hi & DECOMPOSED_IMPLICIT_BIT)) {
         frac_add(a, a, a);
     } else {
         a->exp += 1;
     }
     return;
 
  d_nan:
     float_raise(float_flag_invalid | float_flag_invalid_sqrt, status);
     parts_default_nan(a, status);
 }
 
 /*
  * Rounds the floating-point value `a' to an integer, and returns the
  * result as a floating-point value. The operation is performed
  * according to the IEC/IEEE Standard for Binary Floating-Point
  * Arithmetic.
  *
  * parts_round_to_int_normal is an internal helper function for
  * normal numbers only, returning true for inexact but not directly
  * raising float_flag_inexact.
  */
 static bool partsN(round_to_int_normal)(FloatPartsN *a, FloatRoundMode rmode,
                                         int scale, int frac_size)
 {
     uint64_t frac_lsb, frac_lsbm1, rnd_even_mask, rnd_mask, inc;
     int shift_adj;
 
     scale = MIN(MAX(scale, -0x10000), 0x10000);
     a->exp += scale;
 
     if (a->exp < 0) {
         bool one;
 
         /* All fractional */
         switch (rmode) {
         case float_round_nearest_even:
             one = false;
             if (a->exp == -1) {
                 FloatPartsN tmp;
                 /* Shift left one, discarding DECOMPOSED_IMPLICIT_BIT */
                 frac_add(&tmp, a, a);
                 /* Anything remaining means frac > 0.5. */
                 one = !frac_eqz(&tmp);
             }
             break;
         case float_round_ties_away:
             one = a->exp == -1;
             break;
         case float_round_to_zero:
             one = false;
             break;
         case float_round_up:
             one = !a->sign;
             break;
         case float_round_down:
             one = a->sign;
             break;
         case float_round_to_odd:
             one = true;
             break;
         default:
             g_assert_not_reached();
         }
 
         frac_clear(a);
         a->exp = 0;
         if (one) {
             a->frac_hi = DECOMPOSED_IMPLICIT_BIT;
         } else {
             a->cls = float_class_zero;
         }
         return true;
     }
 
     if (a->exp >= frac_size) {
         /* All integral */
         return false;
     }
 
     if (N > 64 && a->exp < N - 64) {
         /*
          * Rounding is not in the low word -- shift lsb to bit 2,
          * which leaves room for sticky and rounding bit.
          */
         shift_adj = (N - 1) - (a->exp + 2);
         frac_shrjam(a, shift_adj);
         frac_lsb = 1 << 2;
     } else {
         shift_adj = 0;
         frac_lsb = DECOMPOSED_IMPLICIT_BIT >> (a->exp & 63);
     }
 
     frac_lsbm1 = frac_lsb >> 1;
     rnd_mask = frac_lsb - 1;
     rnd_even_mask = rnd_mask | frac_lsb;
 
     if (!(a->frac_lo & rnd_mask)) {
         /* Fractional bits already clear, undo the shift above. */
         frac_shl(a, shift_adj);
         return false;
     }
 
     switch (rmode) {
     case float_round_nearest_even:
         inc = ((a->frac_lo & rnd_even_mask) != frac_lsbm1 ? frac_lsbm1 : 0);
         break;
     case float_round_ties_away:
         inc = frac_lsbm1;
         break;
     case float_round_to_zero:
         inc = 0;
         break;
     case float_round_up:
         inc = a->sign ? 0 : rnd_mask;
         break;
     case float_round_down:
         inc = a->sign ? rnd_mask : 0;
         break;
     case float_round_to_odd:
         inc = a->frac_lo & frac_lsb ? 0 : rnd_mask;
         break;
     default:
         g_assert_not_reached();
     }
 
     if (shift_adj == 0) {
         if (frac_addi(a, a, inc)) {
             frac_shr(a, 1);
             a->frac_hi |= DECOMPOSED_IMPLICIT_BIT;
             a->exp++;
         }
         a->frac_lo &= ~rnd_mask;
     } else {
         frac_addi(a, a, inc);
         a->frac_lo &= ~rnd_mask;
         /* Be careful shifting back, not to overflow */
         frac_shl(a, shift_adj - 1);
         if (a->frac_hi & DECOMPOSED_IMPLICIT_BIT) {
             a->exp++;
         } else {
             frac_add(a, a, a);
         }
     }
     return true;
 }
 
 static void partsN(round_to_int)(FloatPartsN *a, FloatRoundMode rmode,
                                  int scale, float_status *s,
                                  const FloatFmt *fmt)
 {
     switch (a->cls) {
     case float_class_qnan:
     case float_class_snan:
         parts_return_nan(a, s);
         break;
     case float_class_zero:
     case float_class_inf:
         break;
     case float_class_normal:
         if (parts_round_to_int_normal(a, rmode, scale, fmt->frac_size)) {
             float_raise(float_flag_inexact, s);
         }
         break;
     default:
         g_assert_not_reached();
     }
 }
 
 /*
  * Returns the result of converting the floating-point value `a' to
  * the two's complement integer format. The conversion is performed
  * according to the IEC/IEEE Standard for Binary Floating-Point
  * Arithmetic---which means in particular that the conversion is
  * rounded according to the current rounding mode. If `a' is a NaN,
  * the largest positive integer is returned. Otherwise, if the
  * conversion overflows, the largest integer with the same sign as `a'
  * is returned.
  */
 static int64_t partsN(float_to_sint)(FloatPartsN *p, FloatRoundMode rmode,
                                      int scale, int64_t min, int64_t max,
                                      float_status *s)
 {
     int flags = 0;
     uint64_t r;
 
     switch (p->cls) {
     case float_class_snan:
         flags |= float_flag_invalid_snan;
         /* fall through */
     case float_class_qnan:
         flags |= float_flag_invalid;
         r = max;
         break;
 
     case float_class_inf:
         flags = float_flag_invalid | float_flag_invalid_cvti;
         r = p->sign ? min : max;
         break;
 
     case float_class_zero:
         return 0;
 
     case float_class_normal:
         /* TODO: N - 2 is frac_size for rounding; could use input fmt. */
         if (parts_round_to_int_normal(p, rmode, scale, N - 2)) {
             flags = float_flag_inexact;
         }
 
         if (p->exp <= DECOMPOSED_BINARY_POINT) {
             r = p->frac_hi >> (DECOMPOSED_BINARY_POINT - p->exp);
         } else {
             r = UINT64_MAX;
         }
         if (p->sign) {
             if (r <= -(uint64_t)min) {
                 r = -r;
             } else {
                 flags = float_flag_invalid | float_flag_invalid_cvti;
                 r = min;
             }
         } else if (r > max) {
             flags = float_flag_invalid | float_flag_invalid_cvti;
             r = max;
         }
         break;
 
     default:
         g_assert_not_reached();
     }
 
     float_raise(flags, s);
     return r;
 }
 
 /*
  *  Returns the result of converting the floating-point value `a' to
  *  the unsigned integer format. The conversion is performed according
  *  to the IEC/IEEE Standard for Binary Floating-Point
  *  Arithmetic---which means in particular that the conversion is
  *  rounded according to the current rounding mode. If `a' is a NaN,
  *  the largest unsigned integer is returned. Otherwise, if the
  *  conversion overflows, the largest unsigned integer is returned. If
  *  the 'a' is negative, the result is rounded and zero is returned;
  *  values that do not round to zero will raise the inexact exception
  *  flag.
  */
 static uint64_t partsN(float_to_uint)(FloatPartsN *p, FloatRoundMode rmode,
                                       int scale, uint64_t max, float_status *s)
 {
     int flags = 0;
     uint64_t r;
 
     switch (p->cls) {
     case float_class_snan:
         flags |= float_flag_invalid_snan;
         /* fall through */
     case float_class_qnan:
         flags |= float_flag_invalid;
         r = max;
         break;
 
     case float_class_inf:
         flags = float_flag_invalid | float_flag_invalid_cvti;
         r = p->sign ? 0 : max;
         break;
 
     case float_class_zero:
         return 0;
 
     case float_class_normal:
         /* TODO: N - 2 is frac_size for rounding; could use input fmt. */
         if (parts_round_to_int_normal(p, rmode, scale, N - 2)) {
             flags = float_flag_inexact;
             if (p->cls == float_class_zero) {
                 r = 0;
                 break;
             }
         }
 
         if (p->sign) {
             flags = float_flag_invalid | float_flag_invalid_cvti;
             r = 0;
         } else if (p->exp > DECOMPOSED_BINARY_POINT) {
             flags = float_flag_invalid | float_flag_invalid_cvti;
             r = max;
         } else {
             r = p->frac_hi >> (DECOMPOSED_BINARY_POINT - p->exp);
             if (r > max) {
                 flags = float_flag_invalid | float_flag_invalid_cvti;
                 r = max;
             }
         }
         break;
 
     default:
         g_assert_not_reached();
     }
 
     float_raise(flags, s);
     return r;
 }
 
 /*
  * Integer to float conversions
  *
  * Returns the result of converting the two's complement integer `a'
  * to the floating-point format. The conversion is performed according
  * to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
  */
 static void partsN(sint_to_float)(FloatPartsN *p, int64_t a,
                                   int scale, float_status *s)
 {
     uint64_t f = a;
     int shift;
 
     memset(p, 0, sizeof(*p));
 
     if (a == 0) {
         p->cls = float_class_zero;
         return;
     }
 
     p->cls = float_class_normal;
     if (a < 0) {
         f = -f;
         p->sign = true;
     }
     shift = clz64(f);
     scale = MIN(MAX(scale, -0x10000), 0x10000);
 
     p->exp = DECOMPOSED_BINARY_POINT - shift + scale;
     p->frac_hi = f << shift;
 }
 
 /*
  * Unsigned Integer to float conversions
  *
  * Returns the result of converting the unsigned integer `a' to the
  * floating-point format. The conversion is performed according to the
  * IEC/IEEE Standard for Binary Floating-Point Arithmetic.
  */
 static void partsN(uint_to_float)(FloatPartsN *p, uint64_t a,
                                   int scale, float_status *status)
 {
     memset(p, 0, sizeof(*p));
 
     if (a == 0) {
         p->cls = float_class_zero;
     } else {
         int shift = clz64(a);
         scale = MIN(MAX(scale, -0x10000), 0x10000);
         p->cls = float_class_normal;
         p->exp = DECOMPOSED_BINARY_POINT - shift + scale;
         p->frac_hi = a << shift;
     }
 }
 
 /*
  * Float min/max.
  */
 static FloatPartsN *partsN(minmax)(FloatPartsN *a, FloatPartsN *b,
                                    float_status *s, int flags)
 {
     int ab_mask = float_cmask(a->cls) | float_cmask(b->cls);
     int a_exp, b_exp, cmp;
 
     if (unlikely(ab_mask & float_cmask_anynan)) {
         /*
          * For minNum/maxNum (IEEE 754-2008)
          * or minimumNumber/maximumNumber (IEEE 754-2019),
          * if one operand is a QNaN, and the other
          * operand is numerical, then return numerical argument.
          */
         if ((flags & (minmax_isnum | minmax_isnumber))
             && !(ab_mask & float_cmask_snan)
             && (ab_mask & ~float_cmask_qnan)) {
             return is_nan(a->cls) ? b : a;
         }
 
         /*
          * In IEEE 754-2019, minNum, maxNum, minNumMag and maxNumMag
          * are removed and replaced with minimum, minimumNumber, maximum
          * and maximumNumber.
          * minimumNumber/maximumNumber behavior for SNaN is changed to:
          *   If both operands are NaNs, a QNaN is returned.
          *   If either operand is a SNaN,
          *   an invalid operation exception is signaled,
          *   but unless both operands are NaNs,
          *   the SNaN is otherwise ignored and not converted to a QNaN.
          */
         if ((flags & minmax_isnumber)
             && (ab_mask & float_cmask_snan)
             && (ab_mask & ~float_cmask_anynan)) {
             float_raise(float_flag_invalid, s);
             return is_nan(a->cls) ? b : a;
         }
 
         return parts_pick_nan(a, b, s);
     }
 
     a_exp = a->exp;
     b_exp = b->exp;
 
     if (unlikely(ab_mask != float_cmask_normal)) {
         switch (a->cls) {
         case float_class_normal:
             break;
         case float_class_inf:
             a_exp = INT16_MAX;
             break;
         case float_class_zero:
             a_exp = INT16_MIN;
             break;
         default:
             g_assert_not_reached();
             break;
         }
         switch (b->cls) {
         case float_class_normal:
             break;
         case float_class_inf:
             b_exp = INT16_MAX;
             break;
         case float_class_zero:
             b_exp = INT16_MIN;
             break;
         default:
             g_assert_not_reached();
             break;
         }
     }
 
     /* Compare magnitudes. */
     cmp = a_exp - b_exp;
     if (cmp == 0) {
         cmp = frac_cmp(a, b);
     }
 
     /*
      * Take the sign into account.
      * For ismag, only do this if the magnitudes are equal.
      */
     if (!(flags & minmax_ismag) || cmp == 0) {
         if (a->sign != b->sign) {
             /* For differing signs, the negative operand is less. */
             cmp = a->sign ? -1 : 1;
         } else if (a->sign) {
             /* For two negative operands, invert the magnitude comparison. */
             cmp = -cmp;
         }
     }
 
     if (flags & minmax_ismin) {
         cmp = -cmp;
     }
     return cmp < 0 ? b : a;
 }
 
 /*
  * Floating point compare
  */
 static FloatRelation partsN(compare)(FloatPartsN *a, FloatPartsN *b,
                                      float_status *s, bool is_quiet)
 {
     int ab_mask = float_cmask(a->cls) | float_cmask(b->cls);
 
     if (likely(ab_mask == float_cmask_normal)) {
         FloatRelation cmp;
 
         if (a->sign != b->sign) {
             goto a_sign;
         }
         if (a->exp == b->exp) {
             cmp = frac_cmp(a, b);
         } else if (a->exp < b->exp) {
             cmp = float_relation_less;
         } else {
             cmp = float_relation_greater;
         }
         if (a->sign) {
             cmp = -cmp;
         }
         return cmp;
     }
 
     if (unlikely(ab_mask & float_cmask_anynan)) {
         if (ab_mask & float_cmask_snan) {
             float_raise(float_flag_invalid | float_flag_invalid_snan, s);
         } else if (!is_quiet) {
             float_raise(float_flag_invalid, s);
         }
         return float_relation_unordered;
     }
 
     if (ab_mask & float_cmask_zero) {
         if (ab_mask == float_cmask_zero) {
             return float_relation_equal;
         } else if (a->cls == float_class_zero) {
             goto b_sign;
         } else {
             goto a_sign;
         }
     }
 
     if (ab_mask == float_cmask_inf) {
         if (a->sign == b->sign) {
             return float_relation_equal;
         }
     } else if (b->cls == float_class_inf) {
         goto b_sign;
     } else {
         g_assert(a->cls == float_class_inf);
     }
 
  a_sign:
     return a->sign ? float_relation_less : float_relation_greater;
  b_sign:
     return b->sign ? float_relation_greater : float_relation_less;
 }
 
 /*
  * Multiply A by 2 raised to the power N.
  */
 static void partsN(scalbn)(FloatPartsN *a, int n, float_status *s)
 {
     switch (a->cls) {
     case float_class_snan:
     case float_class_qnan:
         parts_return_nan(a, s);
         break;
     case float_class_zero:
     case float_class_inf:
         break;
     case float_class_normal:
         a->exp += MIN(MAX(n, -0x10000), 0x10000);
         break;
     default:
         g_assert_not_reached();
     }
 }
 
 /*
  * Return log2(A)
  */
 static void partsN(log2)(FloatPartsN *a, float_status *s, const FloatFmt *fmt)
 {
     uint64_t a0, a1, r, t, ign;
     FloatPartsN f;
     int i, n, a_exp, f_exp;
 
     if (unlikely(a->cls != float_class_normal)) {
         switch (a->cls) {
         case float_class_snan:
         case float_class_qnan:
             parts_return_nan(a, s);
             return;
         case float_class_zero:
             float_raise(float_flag_divbyzero, s);
             /* log2(0) = -inf */
             a->cls = float_class_inf;
             a->sign = 1;
             return;
         case float_class_inf:
             if (unlikely(a->sign)) {
                 goto d_nan;
             }
             return;
         default:
             break;
         }
         g_assert_not_reached();
     }
     if (unlikely(a->sign)) {
         goto d_nan;
     }
 
     /* TODO: This algorithm looses bits too quickly for float128. */
     g_assert(N == 64);
 
     a_exp = a->exp;
     f_exp = -1;
 
     r = 0;
     t = DECOMPOSED_IMPLICIT_BIT;
     a0 = a->frac_hi;
     a1 = 0;
 
     n = fmt->frac_size + 2;
     if (unlikely(a_exp == -1)) {
         /*
          * When a_exp == -1, we're computing the log2 of a value [0.5,1.0).
          * When the value is very close to 1.0, there are lots of 1's in
          * the msb parts of the fraction.  At the end, when we subtract
          * this value from -1.0, we can see a catastrophic loss of precision,
          * as 0x800..000 - 0x7ff..ffx becomes 0x000..00y, leaving only the
          * bits of y in the final result.  To minimize this, compute as many
          * digits as we can.
          * ??? This case needs another algorithm to avoid this.
          */
         n = fmt->frac_size * 2 + 2;
         /* Don't compute a value overlapping the sticky bit */
         n = MIN(n, 62);
     }
 
     for (i = 0; i < n; i++) {
         if (a1) {
             mul128To256(a0, a1, a0, a1, &a0, &a1, &ign, &ign);
         } else if (a0 & 0xffffffffull) {
             mul64To128(a0, a0, &a0, &a1);
         } else if (a0 & ~DECOMPOSED_IMPLICIT_BIT) {
             a0 >>= 32;
             a0 *= a0;
         } else {
             goto exact;
         }
 
         if (a0 & DECOMPOSED_IMPLICIT_BIT) {
             if (unlikely(a_exp == 0 && r == 0)) {
                 /*
                  * When a_exp == 0, we're computing the log2 of a value
                  * [1.0,2.0).  When the value is very close to 1.0, there
                  * are lots of 0's in the msb parts of the fraction.
                  * We need to compute more digits to produce a correct
                  * result -- restart at the top of the fraction.
                  * ??? This is likely to lose precision quickly, as for
                  * float128; we may need another method.
                  */
                 f_exp -= i;
                 t = r = DECOMPOSED_IMPLICIT_BIT;
                 i = 0;
             } else {
                 r |= t;
             }
         } else {
             add128(a0, a1, a0, a1, &a0, &a1);
         }
         t >>= 1;
     }
 
     /* Set sticky for inexact. */
     r |= (a1 || a0 & ~DECOMPOSED_IMPLICIT_BIT);
 
  exact:
     parts_sint_to_float(a, a_exp, 0, s);
     if (r == 0) {
         return;
     }
 
     memset(&f, 0, sizeof(f));
     f.cls = float_class_normal;
     f.frac_hi = r;
     f.exp = f_exp - frac_normalize(&f);
 
     if (a_exp < 0) {
         parts_sub_normal(a, &f);
     } else if (a_exp > 0) {
         parts_add_normal(a, &f);
     } else {
         *a = f;
     }
     return;
 
  d_nan:
     float_raise(float_flag_invalid, s);
     parts_default_nan(a, s);
 }