Re: [Qemu-devel] [PATCH v3 12/22] fpu/softfloat: re-factor add/sub

[Philippe Mathieu-Daudé]

On 01/24/2018 10:13 AM, Alex Bennée wrote:
> We can now add float16_add/sub and use the common decompose and
> canonicalize functions to have a single implementation for
> float16/32/64 add and sub functions.
> 
> Signed-off-by: Alex Bennée <address@hidden>
> Signed-off-by: Richard Henderson <address@hidden>
> 
> ---
> 
> v3
>  - renames for FloatParts
>  - shorten internal names
>  - fix comment style
>  - use FloatFmt values for for extract/deposit in raw_pack/unpack
>  - move #include bitops from previous patch
>  - add is_nan helper
>  - avoid temp variables where we can
> ---
>  fpu/softfloat.c         | 884 
> +++++++++++++++++++++++++-----------------------
>  include/fpu/softfloat.h |   4 +
>  2 files changed, 461 insertions(+), 427 deletions(-)
> 
> diff --git a/fpu/softfloat.c b/fpu/softfloat.c
> index 568d555595..b179790889 100644
> --- a/fpu/softfloat.c
> +++ b/fpu/softfloat.c
> @@ -83,6 +83,7 @@ this code that are retained.
>   * target-dependent and needs the TARGET_* macros.
>   */
>  #include "qemu/osdep.h"
> +#include "qemu/bitops.h"
>  #include "fpu/softfloat.h"
>  
>  /* We only need stdlib for abort() */
> @@ -270,6 +271,462 @@ static const FloatFmt float64_params = {
>      FLOAT_PARAMS(11, 52)
>  };
>  
> +/* Unpack a float to parts, but do not canonicalize.  */
> +static inline FloatParts unpack_raw(FloatFmt fmt, uint64_t raw)
> +{
> +    const int sign_pos = fmt.frac_size + fmt.exp_size;
> +
> +    return (FloatParts) {
> +        .cls = float_class_unclassified,
> +        .sign = extract64(raw, sign_pos, 1),
> +        .exp = extract64(raw, fmt.frac_size, fmt.exp_size),
> +        .frac = extract64(raw, 0, fmt.frac_size),
> +    };
> +}
> +
> +static inline FloatParts float16_unpack_raw(float16 f)
> +{
> +    return unpack_raw(float16_params, f);
> +}
> +
> +static inline FloatParts float32_unpack_raw(float32 f)
> +{
> +    return unpack_raw(float32_params, f);
> +}
> +
> +static inline FloatParts float64_unpack_raw(float64 f)
> +{
> +    return unpack_raw(float64_params, f);
> +}
> +
> +/* Pack a float from parts, but do not canonicalize.  */
> +static inline uint64_t pack_raw(FloatFmt fmt, FloatParts p)
> +{
> +    const int sign_pos = fmt.frac_size + fmt.exp_size;
> +    uint64_t ret = p.frac;

saving 1 line:

       uint64_t ret = deposit64(p.frac, sign_pos, 1, p.sign);

and then eventually:

       return deposit64(ret, fmt.frac_size, fmt.exp_size, p.exp);

> +
> +    ret = deposit64(ret, fmt.frac_size, fmt.exp_size, p.exp);
> +    ret = deposit64(ret, sign_pos, 1, p.sign);
> +    return ret;
> +}
> +
> +static inline float16 float16_pack_raw(FloatParts p)
> +{
> +    return make_float16(pack_raw(float16_params, p));
> +}
> +
> +static inline float32 float32_pack_raw(FloatParts p)
> +{
> +    return make_float32(pack_raw(float32_params, p));
> +}
> +
> +static inline float64 float64_pack_raw(FloatParts p)
> +{
> +    return make_float64(pack_raw(float64_params, p));
> +}
> +
> +/* Canonicalize EXP and FRAC, setting CLS.  */
> +static FloatParts canonicalize(FloatParts part, const FloatFmt *parm,
> +                               float_status *status)
> +{
> +    if (part.exp == parm->exp_max) {
> +        if (part.frac == 0) {
> +            part.cls = float_class_inf;
> +        } else {
> +#ifdef NO_SIGNALING_NANS
> +            part.cls = float_class_qnan;
> +#else
> +            int64_t msb = part.frac << (parm->frac_shift + 2);
> +            if ((msb < 0) == status->snan_bit_is_one) {
> +                part.cls = float_class_snan;
> +            } else {
> +                part.cls = float_class_qnan;
> +            }
> +#endif
> +        }
> +    } else if (part.exp == 0) {
> +        if (likely(part.frac == 0)) {
> +            part.cls = float_class_zero;
> +        } else if (status->flush_inputs_to_zero) {
> +            float_raise(float_flag_input_denormal, status);
> +            part.cls = float_class_zero;
> +            part.frac = 0;
> +        } else {
> +            int shift = clz64(part.frac) - 1;
> +            part.cls = float_class_normal;
> +            part.exp = parm->frac_shift - parm->exp_bias - shift + 1;
> +            part.frac <<= shift;
> +        }
> +    } else {
> +        part.cls = float_class_normal;
> +        part.exp -= parm->exp_bias;
> +        part.frac = DECOMPOSED_IMPLICIT_BIT + (part.frac << 
> parm->frac_shift);
> +    }
> +    return part;
> +}
> +
> +/* 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 FloatParts round_canonical(FloatParts p, float_status *s,
> +                                  const FloatFmt *parm)
> +{
> +    const uint64_t frac_lsbm1 = parm->frac_lsbm1;
> +    const uint64_t round_mask = parm->round_mask;
> +    const uint64_t roundeven_mask = parm->roundeven_mask;
> +    const int exp_max = parm->exp_max;
> +    const int frac_shift = parm->frac_shift;
> +    uint64_t frac, inc;
> +    int exp, flags = 0;
> +    bool overflow_norm;
> +
> +    frac = p.frac;
> +    exp = p.exp;
> +
> +    switch (p.cls) {
> +    case float_class_normal:
> +        switch (s->float_rounding_mode) {
> +        case float_round_nearest_even:
> +            overflow_norm = false;
> +            inc = ((frac & 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;
> +        default:
> +            g_assert_not_reached();
> +        }
> +
> +        exp += parm->exp_bias;
> +        if (likely(exp > 0)) {
> +            if (frac & round_mask) {
> +                flags |= float_flag_inexact;
> +                frac += inc;
> +                if (frac & DECOMPOSED_OVERFLOW_BIT) {
> +                    frac >>= 1;
> +                    exp++;
> +                }
> +            }
> +            frac >>= frac_shift;
> +
> +            if (unlikely(exp >= exp_max)) {
> +                flags |= float_flag_overflow | float_flag_inexact;
> +                if (overflow_norm) {
> +                    exp = exp_max - 1;
> +                    frac = -1;
> +                } else {
> +                    p.cls = float_class_inf;
> +                    goto do_inf;
> +                }
> +            }
> +        } else if (s->flush_to_zero) {
> +            flags |= float_flag_output_denormal;
> +            p.cls = float_class_zero;
> +            goto do_zero;
> +        } else {
> +            bool is_tiny = (s->float_detect_tininess
> +                            == float_tininess_before_rounding)
> +                        || (exp < 0)
> +                        || !((frac + inc) & DECOMPOSED_OVERFLOW_BIT);
> +
> +            shift64RightJamming(frac, 1 - exp, &frac);
> +            if (frac & round_mask) {
> +                /* Need to recompute round-to-even.  */
> +                if (s->float_rounding_mode == float_round_nearest_even) {
> +                    inc = ((frac & roundeven_mask) != frac_lsbm1
> +                           ? frac_lsbm1 : 0);
> +                }
> +                flags |= float_flag_inexact;
> +                frac += inc;
> +            }
> +
> +            exp = (frac & DECOMPOSED_IMPLICIT_BIT ? 1 : 0);
> +            frac >>= frac_shift;
> +
> +            if (is_tiny && (flags & float_flag_inexact)) {
> +                flags |= float_flag_underflow;
> +            }
> +            if (exp == 0 && frac == 0) {
> +                p.cls = float_class_zero;
> +            }
> +        }
> +        break;
> +
> +    case float_class_zero:
> +    do_zero:
> +        exp = 0;
> +        frac = 0;
> +        break;
> +
> +    case float_class_inf:
> +    do_inf:
> +        exp = exp_max;
> +        frac = 0;
> +        break;
> +
> +    case float_class_qnan:
> +    case float_class_snan:
> +        exp = exp_max;
> +        break;
> +
> +    default:
> +        g_assert_not_reached();
> +    }
> +
> +    float_raise(flags, s);
> +    p.exp = exp;
> +    p.frac = frac;
> +    return p;
> +}
> +
> +static FloatParts float16_unpack_canonical(float16 f, float_status *s)
> +{
> +    return canonicalize(float16_unpack_raw(f), &float16_params, s);
> +}
> +
> +static float16 float16_round_pack_canonical(FloatParts p, float_status *s)
> +{
> +    switch (p.cls) {
> +    case float_class_dnan:
> +        return float16_default_nan(s);
> +    case float_class_msnan:
> +        return float16_maybe_silence_nan(float16_pack_raw(p), s);
> +    default:
> +        p = round_canonical(p, s, &float16_params);
> +        return float16_pack_raw(p);
> +    }
> +}
> +
> +static FloatParts float32_unpack_canonical(float32 f, float_status *s)
> +{
> +    return canonicalize(float32_unpack_raw(f), &float32_params, s);
> +}
> +
> +static float32 float32_round_pack_canonical(FloatParts p, float_status *s)
> +{
> +    switch (p.cls) {
> +    case float_class_dnan:
> +        return float32_default_nan(s);
> +    case float_class_msnan:
> +        return float32_maybe_silence_nan(float32_pack_raw(p), s);
> +    default:
> +        p = round_canonical(p, s, &float32_params);
> +        return float32_pack_raw(p);
> +    }
> +}
> +
> +static FloatParts float64_unpack_canonical(float64 f, float_status *s)
> +{
> +    return canonicalize(float64_unpack_raw(f), &float64_params, s);
> +}
> +
> +static float64 float64_round_pack_canonical(FloatParts p, float_status *s)
> +{
> +    switch (p.cls) {
> +    case float_class_dnan:
> +        return float64_default_nan(s);
> +    case float_class_msnan:
> +        return float64_maybe_silence_nan(float64_pack_raw(p), s);
> +    default:
> +        p = round_canonical(p, s, &float64_params);
> +        return float64_pack_raw(p);
> +    }
> +}
> +
> +static FloatParts pick_nan_parts(FloatParts a, FloatParts b, float_status *s)
> +{
> +    if (a.cls == float_class_snan || b.cls == float_class_snan) {
> +        s->float_exception_flags |= float_flag_invalid;
> +    }
> +
> +    if (s->default_nan_mode) {
> +        a.cls = float_class_dnan;
> +    } else {
> +        if (pickNaN(a.cls == float_class_qnan,
> +                    a.cls == float_class_snan,
> +                    b.cls == float_class_qnan,
> +                    b.cls == float_class_snan,
> +                    a.frac > b.frac
> +                    || (a.frac == b.frac && a.sign < b.sign))) {
> +            a = b;
> +        }
> +        a.cls = float_class_msnan;
> +    }
> +    return a;
> +}
> +
> +/* Simple helper for checking if FloatClass is any NaN */
> +
> +static bool is_nan(FloatClass c)
> +{
> +    return c >= float_class_qnan;
> +}
> +
> +/*
> + * 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 FloatParts addsub_floats(FloatParts a, FloatParts b, bool subtract,
> +                                float_status *s)
> +{
> +    bool a_sign = a.sign;
> +    bool b_sign = b.sign ^ subtract;
> +
> +    if (a_sign != b_sign) {
> +        /* Subtraction */
> +
> +        if (a.cls == float_class_normal && b.cls == float_class_normal) {
> +            if (a.exp > b.exp || (a.exp == b.exp && a.frac >= b.frac)) {
> +                shift64RightJamming(b.frac, a.exp - b.exp, &b.frac);
> +                a.frac = a.frac - b.frac;
> +            } else {
> +                shift64RightJamming(a.frac, b.exp - a.exp, &a.frac);
> +                a.frac = b.frac - a.frac;
> +                a.exp = b.exp;
> +                a_sign ^= 1;
> +            }
> +
> +            if (a.frac == 0) {
> +                a.cls = float_class_zero;
> +                a.sign = s->float_rounding_mode == float_round_down;
> +            } else {
> +                int shift = clz64(a.frac) - 1;
> +                a.frac = a.frac << shift;
> +                a.exp = a.exp - shift;
> +                a.sign = a_sign;
> +            }
> +            return a;
> +        }
> +        if (is_nan(a.cls) || is_nan(b.cls)) {
> +            return pick_nan_parts(a, b, s);
> +        }
> +        if (a.cls == float_class_inf) {
> +            if (b.cls == float_class_inf) {
> +                float_raise(float_flag_invalid, s);
> +                a.cls = float_class_dnan;
> +            }
> +            return a;
> +        }
> +        if (a.cls == float_class_zero && b.cls == float_class_zero) {
> +            a.sign = s->float_rounding_mode == float_round_down;
> +            return a;
> +        }
> +        if (a.cls == float_class_zero || b.cls == float_class_inf) {
> +            b.sign = a_sign ^ 1;
> +            return b;
> +        }
> +        if (b.cls == float_class_zero) {
> +            return a;
> +        }
> +    } else {
> +        /* Addition */
> +        if (a.cls == float_class_normal && b.cls == float_class_normal) {
> +            if (a.exp > b.exp) {
> +                shift64RightJamming(b.frac, a.exp - b.exp, &b.frac);
> +            } else if (a.exp < b.exp) {
> +                shift64RightJamming(a.frac, b.exp - a.exp, &a.frac);
> +                a.exp = b.exp;
> +            }
> +            a.frac += b.frac;
> +            if (a.frac & DECOMPOSED_OVERFLOW_BIT) {
> +                a.frac >>= 1;
> +                a.exp += 1;
> +            }
> +            return a;
> +        }
> +        if (is_nan(a.cls) || is_nan(b.cls)) {
> +            return pick_nan_parts(a, b, s);
> +        }
> +        if (a.cls == float_class_inf || b.cls == float_class_zero) {
> +            return a;
> +        }
> +        if (b.cls == float_class_inf || a.cls == float_class_zero) {
> +            b.sign = b_sign;
> +            return b;
> +        }
> +    }
> +    g_assert_not_reached();
> +}
> +
> +/*
> + * Returns the result of adding or subtracting the floating-point
> + * values `a' and `b'. The operation is performed according to the
> + * IEC/IEEE Standard for Binary Floating-Point Arithmetic.
> + */
> +
> +float16 float16_add(float16 a, float16 b, float_status *status)
> +{
> +    FloatParts pa = float16_unpack_canonical(a, status);
> +    FloatParts pb = float16_unpack_canonical(b, status);
> +    FloatParts pr = addsub_floats(pa, pb, false, status);
> +
> +    return float16_round_pack_canonical(pr, status);
> +}
> +
> +float32 float32_add(float32 a, float32 b, float_status *status)
> +{
> +    FloatParts pa = float32_unpack_canonical(a, status);
> +    FloatParts pb = float32_unpack_canonical(b, status);
> +    FloatParts pr = addsub_floats(pa, pb, false, status);
> +
> +    return float32_round_pack_canonical(pr, status);
> +}
> +
> +float64 float64_add(float64 a, float64 b, float_status *status)
> +{
> +    FloatParts pa = float64_unpack_canonical(a, status);
> +    FloatParts pb = float64_unpack_canonical(b, status);
> +    FloatParts pr = addsub_floats(pa, pb, false, status);
> +
> +    return float64_round_pack_canonical(pr, status);
> +}
> +
> +float16 float16_sub(float16 a, float16 b, float_status *status)
> +{
> +    FloatParts pa = float16_unpack_canonical(a, status);
> +    FloatParts pb = float16_unpack_canonical(b, status);
> +    FloatParts pr = addsub_floats(pa, pb, true, status);
> +
> +    return float16_round_pack_canonical(pr, status);
> +}
> +
> +float32 float32_sub(float32 a, float32 b, float_status *status)
> +{
> +    FloatParts pa = float32_unpack_canonical(a, status);
> +    FloatParts pb = float32_unpack_canonical(b, status);
> +    FloatParts pr = addsub_floats(pa, pb, true, status);
> +
> +    return float32_round_pack_canonical(pr, status);
> +}
> +
> +float64 float64_sub(float64 a, float64 b, float_status *status)
> +{
> +    FloatParts pa = float64_unpack_canonical(a, status);
> +    FloatParts pb = float64_unpack_canonical(b, status);
> +    FloatParts pr = addsub_floats(pa, pb, true, status);
> +
> +    return float64_round_pack_canonical(pr, status);
> +}
> +
>  
> /*----------------------------------------------------------------------------
>  | Takes a 64-bit fixed-point value `absZ' with binary point between bits 6
>  | and 7, and returns the properly rounded 32-bit integer corresponding to the
> @@ -2081,220 +2538,6 @@ float32 float32_round_to_int(float32 a, float_status 
> *status)
>  
>  }
>  
> -/*----------------------------------------------------------------------------
> -| Returns the result of adding the absolute values of the single-precision
> -| floating-point values `a' and `b'.  If `zSign' is 1, the sum is negated
> -| before being returned.  `zSign' is ignored if the result is a NaN.
> -| The addition is performed according to the IEC/IEEE Standard for Binary
> -| Floating-Point Arithmetic.
> -*----------------------------------------------------------------------------*/
> -
> -static float32 addFloat32Sigs(float32 a, float32 b, flag zSign,
> -                              float_status *status)
> -{
> -    int aExp, bExp, zExp;
> -    uint32_t aSig, bSig, zSig;
> -    int expDiff;
> -
> -    aSig = extractFloat32Frac( a );
> -    aExp = extractFloat32Exp( a );
> -    bSig = extractFloat32Frac( b );
> -    bExp = extractFloat32Exp( b );
> -    expDiff = aExp - bExp;
> -    aSig <<= 6;
> -    bSig <<= 6;
> -    if ( 0 < expDiff ) {
> -        if ( aExp == 0xFF ) {
> -            if (aSig) {
> -                return propagateFloat32NaN(a, b, status);
> -            }
> -            return a;
> -        }
> -        if ( bExp == 0 ) {
> -            --expDiff;
> -        }
> -        else {
> -            bSig |= 0x20000000;
> -        }
> -        shift32RightJamming( bSig, expDiff, &bSig );
> -        zExp = aExp;
> -    }
> -    else if ( expDiff < 0 ) {
> -        if ( bExp == 0xFF ) {
> -            if (bSig) {
> -                return propagateFloat32NaN(a, b, status);
> -            }
> -            return packFloat32( zSign, 0xFF, 0 );
> -        }
> -        if ( aExp == 0 ) {
> -            ++expDiff;
> -        }
> -        else {
> -            aSig |= 0x20000000;
> -        }
> -        shift32RightJamming( aSig, - expDiff, &aSig );
> -        zExp = bExp;
> -    }
> -    else {
> -        if ( aExp == 0xFF ) {
> -            if (aSig | bSig) {
> -                return propagateFloat32NaN(a, b, status);
> -            }
> -            return a;
> -        }
> -        if ( aExp == 0 ) {
> -            if (status->flush_to_zero) {
> -                if (aSig | bSig) {
> -                    float_raise(float_flag_output_denormal, status);
> -                }
> -                return packFloat32(zSign, 0, 0);
> -            }
> -            return packFloat32( zSign, 0, ( aSig + bSig )>>6 );
> -        }
> -        zSig = 0x40000000 + aSig + bSig;
> -        zExp = aExp;
> -        goto roundAndPack;
> -    }
> -    aSig |= 0x20000000;
> -    zSig = ( aSig + bSig )<<1;
> -    --zExp;
> -    if ( (int32_t) zSig < 0 ) {
> -        zSig = aSig + bSig;
> -        ++zExp;
> -    }
> - roundAndPack:
> -    return roundAndPackFloat32(zSign, zExp, zSig, status);
> -
> -}
> -
> -/*----------------------------------------------------------------------------
> -| Returns the result of subtracting the absolute values of the single-
> -| precision floating-point values `a' and `b'.  If `zSign' is 1, the
> -| difference is negated before being returned.  `zSign' is ignored if the
> -| result is a NaN.  The subtraction is performed according to the IEC/IEEE
> -| Standard for Binary Floating-Point Arithmetic.
> -*----------------------------------------------------------------------------*/
> -
> -static float32 subFloat32Sigs(float32 a, float32 b, flag zSign,
> -                              float_status *status)
> -{
> -    int aExp, bExp, zExp;
> -    uint32_t aSig, bSig, zSig;
> -    int expDiff;
> -
> -    aSig = extractFloat32Frac( a );
> -    aExp = extractFloat32Exp( a );
> -    bSig = extractFloat32Frac( b );
> -    bExp = extractFloat32Exp( b );
> -    expDiff = aExp - bExp;
> -    aSig <<= 7;
> -    bSig <<= 7;
> -    if ( 0 < expDiff ) goto aExpBigger;
> -    if ( expDiff < 0 ) goto bExpBigger;
> -    if ( aExp == 0xFF ) {
> -        if (aSig | bSig) {
> -            return propagateFloat32NaN(a, b, status);
> -        }
> -        float_raise(float_flag_invalid, status);
> -        return float32_default_nan(status);
> -    }
> -    if ( aExp == 0 ) {
> -        aExp = 1;
> -        bExp = 1;
> -    }
> -    if ( bSig < aSig ) goto aBigger;
> -    if ( aSig < bSig ) goto bBigger;
> -    return packFloat32(status->float_rounding_mode == float_round_down, 0, 
> 0);
> - bExpBigger:
> -    if ( bExp == 0xFF ) {
> -        if (bSig) {
> -            return propagateFloat32NaN(a, b, status);
> -        }
> -        return packFloat32( zSign ^ 1, 0xFF, 0 );
> -    }
> -    if ( aExp == 0 ) {
> -        ++expDiff;
> -    }
> -    else {
> -        aSig |= 0x40000000;
> -    }
> -    shift32RightJamming( aSig, - expDiff, &aSig );
> -    bSig |= 0x40000000;
> - bBigger:
> -    zSig = bSig - aSig;
> -    zExp = bExp;
> -    zSign ^= 1;
> -    goto normalizeRoundAndPack;
> - aExpBigger:
> -    if ( aExp == 0xFF ) {
> -        if (aSig) {
> -            return propagateFloat32NaN(a, b, status);
> -        }
> -        return a;
> -    }
> -    if ( bExp == 0 ) {
> -        --expDiff;
> -    }
> -    else {
> -        bSig |= 0x40000000;
> -    }
> -    shift32RightJamming( bSig, expDiff, &bSig );
> -    aSig |= 0x40000000;
> - aBigger:
> -    zSig = aSig - bSig;
> -    zExp = aExp;
> - normalizeRoundAndPack:
> -    --zExp;
> -    return normalizeRoundAndPackFloat32(zSign, zExp, zSig, status);
> -
> -}
> -
> -/*----------------------------------------------------------------------------
> -| Returns the result of adding the single-precision floating-point values `a'
> -| and `b'.  The operation is performed according to the IEC/IEEE Standard for
> -| Binary Floating-Point Arithmetic.
> -*----------------------------------------------------------------------------*/
> -
> -float32 float32_add(float32 a, float32 b, float_status *status)
> -{
> -    flag aSign, bSign;
> -    a = float32_squash_input_denormal(a, status);
> -    b = float32_squash_input_denormal(b, status);
> -
> -    aSign = extractFloat32Sign( a );
> -    bSign = extractFloat32Sign( b );
> -    if ( aSign == bSign ) {
> -        return addFloat32Sigs(a, b, aSign, status);
> -    }
> -    else {
> -        return subFloat32Sigs(a, b, aSign, status);
> -    }
> -
> -}
> -
> -/*----------------------------------------------------------------------------
> -| Returns the result of subtracting the single-precision floating-point 
> values
> -| `a' and `b'.  The operation is performed according to the IEC/IEEE Standard
> -| for Binary Floating-Point Arithmetic.
> -*----------------------------------------------------------------------------*/
> -
> -float32 float32_sub(float32 a, float32 b, float_status *status)
> -{
> -    flag aSign, bSign;
> -    a = float32_squash_input_denormal(a, status);
> -    b = float32_squash_input_denormal(b, status);
> -
> -    aSign = extractFloat32Sign( a );
> -    bSign = extractFloat32Sign( b );
> -    if ( aSign == bSign ) {
> -        return subFloat32Sigs(a, b, aSign, status);
> -    }
> -    else {
> -        return addFloat32Sigs(a, b, aSign, status);
> -    }
> -
> -}
> -
>  
> /*----------------------------------------------------------------------------
>  | Returns the result of multiplying the single-precision floating-point 
> values
>  | `a' and `b'.  The operation is performed according to the IEC/IEEE Standard
> @@ -3891,219 +4134,6 @@ float64 float64_trunc_to_int(float64 a, float_status 
> *status)
>      return res;
>  }
>  
> -/*----------------------------------------------------------------------------
> -| Returns the result of adding the absolute values of the double-precision
> -| floating-point values `a' and `b'.  If `zSign' is 1, the sum is negated
> -| before being returned.  `zSign' is ignored if the result is a NaN.
> -| The addition is performed according to the IEC/IEEE Standard for Binary
> -| Floating-Point Arithmetic.
> -*----------------------------------------------------------------------------*/
> -
> -static float64 addFloat64Sigs(float64 a, float64 b, flag zSign,
> -                              float_status *status)
> -{
> -    int aExp, bExp, zExp;
> -    uint64_t aSig, bSig, zSig;
> -    int expDiff;
> -
> -    aSig = extractFloat64Frac( a );
> -    aExp = extractFloat64Exp( a );
> -    bSig = extractFloat64Frac( b );
> -    bExp = extractFloat64Exp( b );
> -    expDiff = aExp - bExp;
> -    aSig <<= 9;
> -    bSig <<= 9;
> -    if ( 0 < expDiff ) {
> -        if ( aExp == 0x7FF ) {
> -            if (aSig) {
> -                return propagateFloat64NaN(a, b, status);
> -            }
> -            return a;
> -        }
> -        if ( bExp == 0 ) {
> -            --expDiff;
> -        }
> -        else {
> -            bSig |= LIT64( 0x2000000000000000 );
> -        }
> -        shift64RightJamming( bSig, expDiff, &bSig );
> -        zExp = aExp;
> -    }
> -    else if ( expDiff < 0 ) {
> -        if ( bExp == 0x7FF ) {
> -            if (bSig) {
> -                return propagateFloat64NaN(a, b, status);
> -            }
> -            return packFloat64( zSign, 0x7FF, 0 );
> -        }
> -        if ( aExp == 0 ) {
> -            ++expDiff;
> -        }
> -        else {
> -            aSig |= LIT64( 0x2000000000000000 );
> -        }
> -        shift64RightJamming( aSig, - expDiff, &aSig );
> -        zExp = bExp;
> -    }
> -    else {
> -        if ( aExp == 0x7FF ) {
> -            if (aSig | bSig) {
> -                return propagateFloat64NaN(a, b, status);
> -            }
> -            return a;
> -        }
> -        if ( aExp == 0 ) {
> -            if (status->flush_to_zero) {
> -                if (aSig | bSig) {
> -                    float_raise(float_flag_output_denormal, status);
> -                }
> -                return packFloat64(zSign, 0, 0);
> -            }
> -            return packFloat64( zSign, 0, ( aSig + bSig )>>9 );
> -        }
> -        zSig = LIT64( 0x4000000000000000 ) + aSig + bSig;
> -        zExp = aExp;
> -        goto roundAndPack;
> -    }
> -    aSig |= LIT64( 0x2000000000000000 );
> -    zSig = ( aSig + bSig )<<1;
> -    --zExp;
> -    if ( (int64_t) zSig < 0 ) {
> -        zSig = aSig + bSig;
> -        ++zExp;
> -    }
> - roundAndPack:
> -    return roundAndPackFloat64(zSign, zExp, zSig, status);
> -
> -}
> -
> -/*----------------------------------------------------------------------------
> -| Returns the result of subtracting the absolute values of the double-
> -| precision floating-point values `a' and `b'.  If `zSign' is 1, the
> -| difference is negated before being returned.  `zSign' is ignored if the
> -| result is a NaN.  The subtraction is performed according to the IEC/IEEE
> -| Standard for Binary Floating-Point Arithmetic.
> -*----------------------------------------------------------------------------*/
> -
> -static float64 subFloat64Sigs(float64 a, float64 b, flag zSign,
> -                              float_status *status)
> -{
> -    int aExp, bExp, zExp;
> -    uint64_t aSig, bSig, zSig;
> -    int expDiff;
> -
> -    aSig = extractFloat64Frac( a );
> -    aExp = extractFloat64Exp( a );
> -    bSig = extractFloat64Frac( b );
> -    bExp = extractFloat64Exp( b );
> -    expDiff = aExp - bExp;
> -    aSig <<= 10;
> -    bSig <<= 10;
> -    if ( 0 < expDiff ) goto aExpBigger;
> -    if ( expDiff < 0 ) goto bExpBigger;
> -    if ( aExp == 0x7FF ) {
> -        if (aSig | bSig) {
> -            return propagateFloat64NaN(a, b, status);
> -        }
> -        float_raise(float_flag_invalid, status);
> -        return float64_default_nan(status);
> -    }
> -    if ( aExp == 0 ) {
> -        aExp = 1;
> -        bExp = 1;
> -    }
> -    if ( bSig < aSig ) goto aBigger;
> -    if ( aSig < bSig ) goto bBigger;
> -    return packFloat64(status->float_rounding_mode == float_round_down, 0, 
> 0);
> - bExpBigger:
> -    if ( bExp == 0x7FF ) {
> -        if (bSig) {
> -            return propagateFloat64NaN(a, b, status);
> -        }
> -        return packFloat64( zSign ^ 1, 0x7FF, 0 );
> -    }
> -    if ( aExp == 0 ) {
> -        ++expDiff;
> -    }
> -    else {
> -        aSig |= LIT64( 0x4000000000000000 );
> -    }
> -    shift64RightJamming( aSig, - expDiff, &aSig );
> -    bSig |= LIT64( 0x4000000000000000 );
> - bBigger:
> -    zSig = bSig - aSig;
> -    zExp = bExp;
> -    zSign ^= 1;
> -    goto normalizeRoundAndPack;
> - aExpBigger:
> -    if ( aExp == 0x7FF ) {
> -        if (aSig) {
> -            return propagateFloat64NaN(a, b, status);
> -        }
> -        return a;
> -    }
> -    if ( bExp == 0 ) {
> -        --expDiff;
> -    }
> -    else {
> -        bSig |= LIT64( 0x4000000000000000 );
> -    }
> -    shift64RightJamming( bSig, expDiff, &bSig );
> -    aSig |= LIT64( 0x4000000000000000 );
> - aBigger:
> -    zSig = aSig - bSig;
> -    zExp = aExp;
> - normalizeRoundAndPack:
> -    --zExp;
> -    return normalizeRoundAndPackFloat64(zSign, zExp, zSig, status);
> -
> -}
> -
> -/*----------------------------------------------------------------------------
> -| Returns the result of adding the double-precision floating-point values `a'
> -| and `b'.  The operation is performed according to the IEC/IEEE Standard for
> -| Binary Floating-Point Arithmetic.
> -*----------------------------------------------------------------------------*/
> -
> -float64 float64_add(float64 a, float64 b, float_status *status)
> -{
> -    flag aSign, bSign;
> -    a = float64_squash_input_denormal(a, status);
> -    b = float64_squash_input_denormal(b, status);
> -
> -    aSign = extractFloat64Sign( a );
> -    bSign = extractFloat64Sign( b );
> -    if ( aSign == bSign ) {
> -        return addFloat64Sigs(a, b, aSign, status);
> -    }
> -    else {
> -        return subFloat64Sigs(a, b, aSign, status);
> -    }
> -
> -}
> -
> -/*----------------------------------------------------------------------------
> -| Returns the result of subtracting the double-precision floating-point 
> values
> -| `a' and `b'.  The operation is performed according to the IEC/IEEE Standard
> -| for Binary Floating-Point Arithmetic.
> -*----------------------------------------------------------------------------*/
> -
> -float64 float64_sub(float64 a, float64 b, float_status *status)
> -{
> -    flag aSign, bSign;
> -    a = float64_squash_input_denormal(a, status);
> -    b = float64_squash_input_denormal(b, status);
> -
> -    aSign = extractFloat64Sign( a );
> -    bSign = extractFloat64Sign( b );
> -    if ( aSign == bSign ) {
> -        return subFloat64Sigs(a, b, aSign, status);
> -    }
> -    else {
> -        return addFloat64Sigs(a, b, aSign, status);
> -    }
> -
> -}
>  
>  
> /*----------------------------------------------------------------------------
>  | Returns the result of multiplying the double-precision floating-point 
> values
> diff --git a/include/fpu/softfloat.h b/include/fpu/softfloat.h
> index 23824a3000..693ece0974 100644
> --- a/include/fpu/softfloat.h
> +++ b/include/fpu/softfloat.h
> @@ -236,6 +236,10 @@ float64 float16_to_float64(float16 a, flag ieee, 
> float_status *status);
>  
> /*----------------------------------------------------------------------------
>  | Software half-precision operations.
>  
> *----------------------------------------------------------------------------*/
> +
> +float16 float16_add(float16, float16, float_status *status);
> +float16 float16_sub(float16, float16, float_status *status);
> +
>  int float16_is_quiet_nan(float16, float_status *status);
>  int float16_is_signaling_nan(float16, float_status *status);
>  float16 float16_maybe_silence_nan(float16, float_status *status);
> 

(note for myself, since I need more time for those 3)
all bug round_canonical() pick_nan_parts() addsub_floats():
Reviewed-by: Philippe Mathieu-Daudé <address@hidden>