diff options
Diffstat (limited to 'fpu/softfloat-parts.c.inc')
| -rw-r--r-- | fpu/softfloat-parts.c.inc | 817 |
1 files changed, 817 insertions, 0 deletions
diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc new file mode 100644 index 0000000000..a897a5a743 --- /dev/null +++ b/fpu/softfloat-parts.c.inc @@ -0,0 +1,817 @@ +/* + * 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, 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, 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, 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)(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 frac_lsb = fmt->frac_lsb; + const uint64_t frac_lsbm1 = fmt->frac_lsbm1; + const uint64_t round_mask = fmt->round_mask; + const uint64_t roundeven_mask = fmt->roundeven_mask; + uint64_t inc; + bool overflow_norm; + int exp, flags = 0; + + if (unlikely(p->cls != float_class_normal)) { + 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(); + } + + switch (s->float_rounding_mode) { + case float_round_nearest_even: + overflow_norm = false; + inc = ((p->frac_lo & roundeven_mask) != frac_lsbm1 ? frac_lsbm1 : 0); + break; + case float_round_ties_away: + overflow_norm = false; + 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; + 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++; + } + } + frac_shr(p, frac_shift); + + 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); + } + } else if (unlikely(exp >= exp_max)) { + flags |= float_flag_overflow | float_flag_inexact; + if (overflow_norm) { + exp = exp_max - 1; + frac_allones(p); + } else { + p->cls = float_class_inf; + exp = exp_max; + frac_clear(p); + } + } + } 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: + inc = ((p->frac_lo & roundeven_mask) != frac_lsbm1 + ? frac_lsbm1 : 0); + break; + case float_round_to_odd: + inc = p->frac_lo & frac_lsb ? 0 : round_mask; + break; + default: + break; + } + flags |= float_flag_inexact; + frac_addi(p, p, inc); + } + + 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); +} + +/* + * 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, 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, 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)) { + goto d_nan; + } + + if (ab_mask & float_cmask_inf) { + if (c->cls == float_class_inf && a->sign != c->sign) { + 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: + float_raise(float_flag_invalid, s); + 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) || + unlikely(ab_mask == float_cmask_inf)) { + float_raise(float_flag_invalid, s); + parts_default_nan(a, s); + return a; + } + + /* 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; +} + +/* + * 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: + case float_class_qnan: + flags = float_flag_invalid; + r = max; + break; + + case float_class_inf: + flags = float_flag_invalid; + 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; + r = min; + } + } else if (r > max) { + flags = float_flag_invalid; + r = max; + } + break; + + default: + g_assert_not_reached(); + } + + float_raise(flags, s); + return r; +} |