diff mbox series

[RFC,15/30] softfloat: half-precision add/sub/mul/div support

Message ID 20171013162438.32458-16-alex.bennee@linaro.org
State New
Headers show
Series v8.2 half-precision support (work-in-progress) | expand

Commit Message

Alex Bennée Oct. 13, 2017, 4:24 p.m. UTC 
Rather than following the SoftFloat3 implementation I've used the same
basic template as the rest of our softfloat code. One minor difference
is the 32bit intermediates end up with the binary point in the same
place as the 32 bit version so the change isn't totally mechanical.

Signed-off-by: Alex Bennée <alex.bennee@linaro.org>
---
 fpu/softfloat.c         | 352 ++++++++++++++++++++++++++++++++++++++++++++++++
 include/fpu/softfloat.h |   6 +
 2 files changed, 358 insertions(+)

Comments

Richard Henderson Oct. 16, 2017, 10:01 p.m. UTC | #1 
On 10/13/2017 09:24 AM, Alex Bennée wrote:
> +static float16 addFloat16Sigs(float16 a, float16 b, flag zSign,
> +                              float_status *status)
> +{
> +    int aExp, bExp, zExp;
> +    uint16_t aSig, bSig, zSig;
> +    int expDiff;
> +
> +    aSig = extractFloat16Frac( a );
> +    aExp = extractFloat16Exp( a );
> +    bSig = extractFloat16Frac( b );
> +    bExp = extractFloat16Exp( b );
> +    expDiff = aExp - bExp;
> +    aSig <<= 3;
> +    bSig <<= 3;
> +    if ( 0 < expDiff ) {
> +        if ( aExp == 0x1F ) {
> +            if (aSig) {
> +                return propagateFloat16NaN(a, b, status);
> +            }
> +            return a;
> +        }
> +        if ( bExp == 0 ) {
> +            --expDiff;
> +        }
> +        else {
> +            bSig |= 0x20000000;
> +        }

This number is wrong.  The 32-bit number is the implicit bit (0x800000) << 6
(which is the shift applied there).  Here it should be 0x400 << 3.  And it
probably should be written that way for clarity.

> +    else if ( expDiff < 0 ) {
> +        if ( bExp == 0x1F ) {
> +            if (bSig) {
> +                return propagateFloat16NaN(a, b, status);
> +            }
> +            return packFloat16( zSign, 0x1F, 0 );
> +        }
> +        if ( aExp == 0 ) {
> +            ++expDiff;
> +        }
> +        else {
> +            aSig |= 0x0400;

This number is wrong because you didn't apply the << 3.

> +        if ( aExp == 0 ) {
> +            if (status->flush_to_zero) {
> +                if (aSig | bSig) {
> +                    float_raise(float_flag_output_denormal, status);
> +                }
> +                return packFloat16(zSign, 0, 0);
> +            }
> +            return packFloat16( zSign, 0, ( aSig + bSig )>>3 );
> +        }
> +        zSig = 0x0400 + aSig + bSig;
> +        zExp = aExp;
> +        goto roundAndPack;
> +    }
> +    aSig |= 0x0400;

Needs << 3.

> +    zSig = ( aSig + bSig )<<1;
> +    --zExp;
> +    if ( (int16_t) zSig < 0 ) {
> +        zSig = aSig + bSig;
> +        ++zExp;
> +    }

As a side-note, I really wish the existing code didn't try to half-optimize
this.  Both aSig and bSig really should have their implicit bits added.  Of
course, that affects this weird zExp dance as well.

> +    int aExp, bExp, zExp;
> +    uint16_t aSig, bSig, zSig;
> +    int expDiff;
> +
> +    aSig = extractFloat16Frac( a );
> +    aExp = extractFloat16Exp( a );
> +    bSig = extractFloat16Frac( b );
> +    bExp = extractFloat16Exp( b );
> +    expDiff = aExp - bExp;
> +    aSig <<= 7;
> +    bSig <<= 7;

You've just lost one frac bit off the left of the uint16_t.

You cannot add that many bias bits and still use that as an intermediate data
type.  You can only shift by 3, as you did for add.  Or, widen to uint32_t.

> +    if ( 0 < expDiff ) goto aExpBigger;
> +    if ( expDiff < 0 ) goto bExpBigger;
> +    if ( aExp == 0xFF ) {
> +        if (aSig | bSig) {
> +            return propagateFloat16NaN(a, b, status);
> +        }
> +        float_raise(float_flag_invalid, status);
> +        return float16_default_nan(status);
> +    }
> +    if ( aExp == 0 ) {
> +        aExp = 1;
> +        bExp = 1;
> +    }
> +    if ( bSig < aSig ) goto aBigger;
> +    if ( aSig < bSig ) goto bBigger;
> +    return packFloat16(status->float_rounding_mode == float_round_down, 0, 0);
> + bExpBigger:
> +    if ( bExp == 0xFF ) {
> +        if (bSig) {
> +            return propagateFloat16NaN(a, b, status);
> +        }
> +        return packFloat16( zSign ^ 1, 0xFF, 0 );
> +    }
> +    if ( aExp == 0 ) {
> +        ++expDiff;
> +    }
> +    else {
> +        aSig |= 0x40000000;
> +    }
> +    shift16RightJamming( aSig, - expDiff, &aSig );
> +    bSig |= 0x40000000;

These should be 0x400 << shl.  It should be obvious that 0x40000000 doesn't fit
in uint16_t.

> + bBigger:
> +    zSig = bSig - aSig;
> +    zExp = bExp;
> +    zSign ^= 1;
> +    goto normalizeRoundAndPack;
> + aExpBigger:
> +    if ( aExp == 0xFF ) {
> +        if (aSig) {
> +            return propagateFloat16NaN(a, b, status);
> +        }
> +        return a;
> +    }
> +    if ( bExp == 0 ) {
> +        --expDiff;
> +    }
> +    else {
> +        bSig |= 0x40000000;
> +    }
> +    shift16RightJamming( bSig, expDiff, &bSig );
> +    aSig |= 0x40000000;

Likewise.

> + aBigger:
> +    zSig = aSig - bSig;
> +    zExp = aExp;
> + normalizeRoundAndPack:
> +    --zExp;
> +    return normalizeRoundAndPackFloat16(zSign, zExp, zSig, status);

This reinforces my point about not using uint16_t

> +float16 float16_mul(float16 a, float16 b, float_status *status)
...
> +    if ( aExp == 0 ) {
> +        if ( aSig == 0 ) return packFloat16( zSign, 0, 0 );
> +        normalizeFloat16Subnormal( aSig, &aExp, &aSig );
> +    }
> +    if ( bExp == 0 ) {
> +        if ( bSig == 0 ) return packFloat16( zSign, 0, 0 );
> +        normalizeFloat16Subnormal( bSig, &bExp, &bSig );
> +    }
> +    zExp = aExp + bExp - 0xF;
> +    /* Add implicit bit */
> +    aSig = ( aSig | 0x0400 )<<4;
> +    bSig = ( bSig | 0x0400 )<<5;

The implicit bit shouldn't be added for denormals.  It *probably* doesn't
matter since the normalization step should be setting the same bit, but it's
certainly less than clear.

Is there a particular reason to add the 9 bits of shifting?  I thought we
demonstrated to your satisfaction that it wasn't required.

> +    /* Max (format " => 0x%x" (* (lsh #x400 4)  (lsh #x400 5))) => 0x20000000
> +     * So shift so binary point from 30/29 to 23/22
> +     */

That's the minimum, not maximum.  Maximum would be

  0x7ff<<4 * 0x7ff<<5 => 0x7fe00200

I am totally unclear why the binary point should be at 23/22.  That doesn't
seem to match up with the shifts you've used for either addition or
subtraction.  Both of which were different, btw...


> +    shift32RightJamming( ( (uint32_t) aSig ) * bSig, 7, &zSig32 );
> +    /* At this point the significand is at the same point as
> +     * float32_mul, so we can do the same test */
> +    if ( 0 <= (int32_t) ( zSig32<<1 ) ) {
> +        zSig32 <<= 1;
> +        --zExp;
> +    }

These two are definitely in the wrong order.  You can certainly normalize
0x20000000 to 0x40000000 via that IF, but not if you've already right-shifted
by 7 bits.


> +/*----------------------------------------------------------------------------
> +| Returns the result of dividing the half-precision floating-point value `a'
> +| by the corresponding value `b'.  The operation is performed according to the
> +| IEC/IEEE Standard for Binary Floating-Point Arithmetic.
> +*----------------------------------------------------------------------------*/
> +
> +float16 float16_div(float16 a, float16 b, float_status *status)
> +{
> +    flag aSign, bSign, zSign;
> +    int aExp, bExp, zExp;
> +    uint32_t aSig, bSig, zSig;
> +    a = float16_squash_input_denormal(a, status);
> +    b = float16_squash_input_denormal(b, status);
> +
> +    aSig = extractFloat16Frac( a );
> +    aExp = extractFloat16Exp( a );
> +    aSign = extractFloat16Sign( a );
> +    bSig = extractFloat16Frac( b );
> +    bExp = extractFloat16Exp( b );
> +    bSign = extractFloat16Sign( b );
> +    zSign = aSign ^ bSign;
> +    if ( aExp == 0xFF ) {
> +        if (aSig) {
> +            return propagateFloat16NaN(a, b, status);
> +        }
> +        if ( bExp == 0xFF ) {
> +            if (bSig) {
> +                return propagateFloat16NaN(a, b, status);
> +            }
> +            float_raise(float_flag_invalid, status);
> +            return float16_default_nan(status);
> +        }
> +        return packFloat16( zSign, 0xFF, 0 );

All these 0xff should be 0x1f for float16.

> +    if ( bExp == 0 ) {
> +        if ( bSig == 0 ) {
> +            if ( ( aExp | aSig ) == 0 ) {
> +                float_raise(float_flag_invalid, status);
> +                return float16_default_nan(status);
> +            }
> +            float_raise(float_flag_divbyzero, status);
> +            return packFloat16( zSign, 0xFF, 0 );
> +        }
> +        normalizeFloat16Subnormal( bSig, &bExp, &bSig );
> +    }
> +    if ( aExp == 0 ) {
> +        if ( aSig == 0 ) return packFloat16( zSign, 0, 0 );
> +        normalizeFloat16Subnormal( aSig, &aExp, &aSig );
> +    }
> +    zExp = aExp - bExp + 0x7D;

The exponent bias correction is wrong.

> +    aSig = ( aSig | 0x00800000 )<<7;
> +    bSig = ( bSig | 0x00800000 )<<8;

Again, the implicit bit isn't correct, and probably should be explicitly set
for subnormals.

> +    if ( bSig <= ( aSig + aSig ) ) {
> +        aSig >>= 1;
> +        ++zExp;
> +    }
> +    zSig = ( ( (uint64_t) aSig )<<16 ) / bSig;

If you hadn't already shifted by 7 & 8, you wouldn't need to make this a
uint64_t for division.  Indeed, those shift were chosen for float32 so that the
float32 fraction was at the top of the uint32_t.  We then normalize the
operands such that result must be fractional.

> +    if ( ( zSig & 0x3F ) == 0 ) {
> +        zSig |= ( (uint64_t) bSig * zSig != ( (uint64_t) aSig )<<16 );
> +    }

Well, the 0x3f isn't right.  I think ideally we'd use div and mod
unconditionally and let the compiler dtrt, e.g.

  aSig <<= 16;
  q = aSig / bSig;
  r = aSig % bSig;
  zSig = q | (r != 0);


r~
diff mbox series

Patch

diff --git a/fpu/softfloat.c b/fpu/softfloat.c
index cf7bf6d4f4..ff967f5525 100644
--- a/fpu/softfloat.c
+++ b/fpu/softfloat.c
@@ -3532,6 +3532,358 @@  static void normalizeFloat16Subnormal(uint32_t aSig, int *zExpPtr,
     *zExpPtr = 1 - shiftCount;
 }
 
+/*----------------------------------------------------------------------------
+| Returns the result of adding the absolute values of the half-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 float16 addFloat16Sigs(float16 a, float16 b, flag zSign,
+                              float_status *status)
+{
+    int aExp, bExp, zExp;
+    uint16_t aSig, bSig, zSig;
+    int expDiff;
+
+    aSig = extractFloat16Frac( a );
+    aExp = extractFloat16Exp( a );
+    bSig = extractFloat16Frac( b );
+    bExp = extractFloat16Exp( b );
+    expDiff = aExp - bExp;
+    aSig <<= 3;
+    bSig <<= 3;
+    if ( 0 < expDiff ) {
+        if ( aExp == 0x1F ) {
+            if (aSig) {
+                return propagateFloat16NaN(a, b, status);
+            }
+            return a;
+        }
+        if ( bExp == 0 ) {
+            --expDiff;
+        }
+        else {
+            bSig |= 0x20000000;
+        }
+        shift16RightJamming( bSig, expDiff, &bSig );
+        zExp = aExp;
+    }
+    else if ( expDiff < 0 ) {
+        if ( bExp == 0x1F ) {
+            if (bSig) {
+                return propagateFloat16NaN(a, b, status);
+            }
+            return packFloat16( zSign, 0x1F, 0 );
+        }
+        if ( aExp == 0 ) {
+            ++expDiff;
+        }
+        else {
+            aSig |= 0x0400;
+        }
+        shift16RightJamming( aSig, - expDiff, &aSig );
+        zExp = bExp;
+    }
+    else {
+        if ( aExp == 0x1F ) {
+            if (aSig | bSig) {
+                return propagateFloat16NaN(a, b, status);
+            }
+            return a;
+        }
+        if ( aExp == 0 ) {
+            if (status->flush_to_zero) {
+                if (aSig | bSig) {
+                    float_raise(float_flag_output_denormal, status);
+                }
+                return packFloat16(zSign, 0, 0);
+            }
+            return packFloat16( zSign, 0, ( aSig + bSig )>>3 );
+        }
+        zSig = 0x0400 + aSig + bSig;
+        zExp = aExp;
+        goto roundAndPack;
+    }
+    aSig |= 0x0400;
+    zSig = ( aSig + bSig )<<1;
+    --zExp;
+    if ( (int16_t) zSig < 0 ) {
+        zSig = aSig + bSig;
+        ++zExp;
+    }
+ roundAndPack:
+    return roundAndPackFloat16(zSign, zExp, zSig, true, status);
+
+}
+
+/*----------------------------------------------------------------------------
+| Returns the result of subtracting the absolute values of the half-
+| 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 float16 subFloat16Sigs(float16 a, float16 b, flag zSign,
+                              float_status *status)
+{
+    int aExp, bExp, zExp;
+    uint16_t aSig, bSig, zSig;
+    int expDiff;
+
+    aSig = extractFloat16Frac( a );
+    aExp = extractFloat16Exp( a );
+    bSig = extractFloat16Frac( b );
+    bExp = extractFloat16Exp( 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 propagateFloat16NaN(a, b, status);
+        }
+        float_raise(float_flag_invalid, status);
+        return float16_default_nan(status);
+    }
+    if ( aExp == 0 ) {
+        aExp = 1;
+        bExp = 1;
+    }
+    if ( bSig < aSig ) goto aBigger;
+    if ( aSig < bSig ) goto bBigger;
+    return packFloat16(status->float_rounding_mode == float_round_down, 0, 0);
+ bExpBigger:
+    if ( bExp == 0xFF ) {
+        if (bSig) {
+            return propagateFloat16NaN(a, b, status);
+        }
+        return packFloat16( zSign ^ 1, 0xFF, 0 );
+    }
+    if ( aExp == 0 ) {
+        ++expDiff;
+    }
+    else {
+        aSig |= 0x40000000;
+    }
+    shift16RightJamming( aSig, - expDiff, &aSig );
+    bSig |= 0x40000000;
+ bBigger:
+    zSig = bSig - aSig;
+    zExp = bExp;
+    zSign ^= 1;
+    goto normalizeRoundAndPack;
+ aExpBigger:
+    if ( aExp == 0xFF ) {
+        if (aSig) {
+            return propagateFloat16NaN(a, b, status);
+        }
+        return a;
+    }
+    if ( bExp == 0 ) {
+        --expDiff;
+    }
+    else {
+        bSig |= 0x40000000;
+    }
+    shift16RightJamming( bSig, expDiff, &bSig );
+    aSig |= 0x40000000;
+ aBigger:
+    zSig = aSig - bSig;
+    zExp = aExp;
+ normalizeRoundAndPack:
+    --zExp;
+    return normalizeRoundAndPackFloat16(zSign, zExp, zSig, status);
+
+}
+
+/*----------------------------------------------------------------------------
+| Returns the result of adding the half-precision 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)
+{
+    flag aSign, bSign;
+    a = float16_squash_input_denormal(a, status);
+    b = float16_squash_input_denormal(b, status);
+
+    aSign = extractFloat16Sign( a );
+    bSign = extractFloat16Sign( b );
+    if ( aSign == bSign ) {
+        return addFloat16Sigs(a, b, aSign, status);
+    }
+    else {
+        return subFloat16Sigs(a, b, aSign, status);
+    }
+
+}
+
+/*----------------------------------------------------------------------------
+| Returns the result of subtracting the half-precision floating-point values
+| `a' and `b'.  The operation is performed according to the IEC/IEEE Standard
+| for Binary Floating-Point Arithmetic.
+*----------------------------------------------------------------------------*/
+
+float16 float16_sub(float16 a, float16 b, float_status *status)
+{
+    flag aSign, bSign;
+    a = float16_squash_input_denormal(a, status);
+    b = float16_squash_input_denormal(b, status);
+
+    aSign = extractFloat16Sign( a );
+    bSign = extractFloat16Sign( b );
+    if ( aSign == bSign ) {
+        return subFloat16Sigs(a, b, aSign, status);
+    }
+    else {
+        return addFloat16Sigs(a, b, aSign, status);
+    }
+
+}
+
+/*----------------------------------------------------------------------------
+| Returns the result of multiplying the half-precision floating-point values
+| `a' and `b'.  The operation is performed according to the IEC/IEEE Standard
+| for Binary Floating-Point Arithmetic.
+*----------------------------------------------------------------------------*/
+
+float16 float16_mul(float16 a, float16 b, float_status *status)
+{
+    flag aSign, bSign, zSign;
+    int aExp, bExp, zExp;
+    uint32_t aSig, bSig;
+    uint32_t zSig32; /* no zSig as zSig32 passed into rp&f */
+
+    a = float16_squash_input_denormal(a, status);
+    b = float16_squash_input_denormal(b, status);
+
+    aSig = extractFloat16Frac( a );
+    aExp = extractFloat16Exp( a );
+    aSign = extractFloat16Sign( a );
+    bSig = extractFloat16Frac( b );
+    bExp = extractFloat16Exp( b );
+    bSign = extractFloat16Sign( b );
+    zSign = aSign ^ bSign;
+    if ( aExp == 0x1F ) {
+        if ( aSig || ( ( bExp == 0x1F ) && bSig ) ) {
+            return propagateFloat16NaN(a, b, status);
+        }
+        if ( ( bExp | bSig ) == 0 ) {
+            float_raise(float_flag_invalid, status);
+            return float16_default_nan(status);
+        }
+        return packFloat16( zSign, 0x1F, 0 );
+    }
+    if ( bExp == 0x1F ) {
+        if (bSig) {
+            return propagateFloat16NaN(a, b, status);
+        }
+        if ( ( aExp | aSig ) == 0 ) {
+            float_raise(float_flag_invalid, status);
+            return float16_default_nan(status);
+        }
+        return packFloat16( zSign, 0x1F, 0 );
+    }
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return packFloat16( zSign, 0, 0 );
+        normalizeFloat16Subnormal( aSig, &aExp, &aSig );
+    }
+    if ( bExp == 0 ) {
+        if ( bSig == 0 ) return packFloat16( zSign, 0, 0 );
+        normalizeFloat16Subnormal( bSig, &bExp, &bSig );
+    }
+    zExp = aExp + bExp - 0xF;
+    /* Add implicit bit */
+    aSig = ( aSig | 0x0400 )<<4;
+    bSig = ( bSig | 0x0400 )<<5;
+    /* Max (format " => 0x%x" (* (lsh #x400 4)  (lsh #x400 5))) => 0x20000000
+     * So shift so binary point from 30/29 to 23/22
+     */
+    shift32RightJamming( ( (uint32_t) aSig ) * bSig, 7, &zSig32 );
+    /* At this point the significand is at the same point as
+     * float32_mul, so we can do the same test */
+    if ( 0 <= (int32_t) ( zSig32<<1 ) ) {
+        zSig32 <<= 1;
+        --zExp;
+    }
+    return roundAndPackFloat16(zSign, zExp, zSig32, true, status);
+}
+
+/*----------------------------------------------------------------------------
+| Returns the result of dividing the half-precision floating-point value `a'
+| by the corresponding value `b'.  The operation is performed according to the
+| IEC/IEEE Standard for Binary Floating-Point Arithmetic.
+*----------------------------------------------------------------------------*/
+
+float16 float16_div(float16 a, float16 b, float_status *status)
+{
+    flag aSign, bSign, zSign;
+    int aExp, bExp, zExp;
+    uint32_t aSig, bSig, zSig;
+    a = float16_squash_input_denormal(a, status);
+    b = float16_squash_input_denormal(b, status);
+
+    aSig = extractFloat16Frac( a );
+    aExp = extractFloat16Exp( a );
+    aSign = extractFloat16Sign( a );
+    bSig = extractFloat16Frac( b );
+    bExp = extractFloat16Exp( b );
+    bSign = extractFloat16Sign( b );
+    zSign = aSign ^ bSign;
+    if ( aExp == 0xFF ) {
+        if (aSig) {
+            return propagateFloat16NaN(a, b, status);
+        }
+        if ( bExp == 0xFF ) {
+            if (bSig) {
+                return propagateFloat16NaN(a, b, status);
+            }
+            float_raise(float_flag_invalid, status);
+            return float16_default_nan(status);
+        }
+        return packFloat16( zSign, 0xFF, 0 );
+    }
+    if ( bExp == 0xFF ) {
+        if (bSig) {
+            return propagateFloat16NaN(a, b, status);
+        }
+        return packFloat16( zSign, 0, 0 );
+    }
+    if ( bExp == 0 ) {
+        if ( bSig == 0 ) {
+            if ( ( aExp | aSig ) == 0 ) {
+                float_raise(float_flag_invalid, status);
+                return float16_default_nan(status);
+            }
+            float_raise(float_flag_divbyzero, status);
+            return packFloat16( zSign, 0xFF, 0 );
+        }
+        normalizeFloat16Subnormal( bSig, &bExp, &bSig );
+    }
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return packFloat16( zSign, 0, 0 );
+        normalizeFloat16Subnormal( aSig, &aExp, &aSig );
+    }
+    zExp = aExp - bExp + 0x7D;
+    aSig = ( aSig | 0x00800000 )<<7;
+    bSig = ( bSig | 0x00800000 )<<8;
+    if ( bSig <= ( aSig + aSig ) ) {
+        aSig >>= 1;
+        ++zExp;
+    }
+    zSig = ( ( (uint64_t) aSig )<<16 ) / bSig;
+    if ( ( zSig & 0x3F ) == 0 ) {
+        zSig |= ( (uint64_t) bSig * zSig != ( (uint64_t) aSig )<<16 );
+    }
+    return roundAndPackFloat16(zSign, zExp, zSig, true, status);
+
+}
+
 /* Half precision floats come in two formats: standard IEEE and "ARM" format.
    The latter gains extra exponent range by omitting the NaN/Inf encodings.  */
 
diff --git a/include/fpu/softfloat.h b/include/fpu/softfloat.h
index d89fdf7675..f1d79b6d03 100644
--- a/include/fpu/softfloat.h
+++ b/include/fpu/softfloat.h
@@ -345,6 +345,12 @@  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);
+float16 float16_mul(float16, float16, float_status *status);
+float16 float16_div(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);