This is one of the simpler manipulations you could make to a floating
point number.

Signed-off-by: Alex Bennée <alex.bennee@linaro.org>
Reviewed-by: Richard Henderson <richard.henderson@linaro.org>

---
v3
  - fix renames
---
 fpu/softfloat.c         | 104 +++++++++++++++---------------------------------
 include/fpu/softfloat.h |   1 +
 2 files changed, 32 insertions(+), 73 deletions(-)

diff --git a/fpu/softfloat.c b/fpu/softfloat.c
index 7baeb0e22d..24b7de30ce 100644
--- a/fpu/softfloat.c
+++ b/fpu/softfloat.c
@@ -1646,6 +1646,37 @@ float64 uint16_to_float64(uint16_t a, float_status *status)
     return uint64_to_float64(a, status);
 }
 
+/* Multiply A by 2 raised to the power N.  */
+static FloatParts scalbn_decomposed(FloatParts a, int n,
+                                          float_status *s)
+{
+    if (a.cls == float_class_normal) {
+        a.exp += n;
+    }
+    return a;
+}
+
+float16 float16_scalbn(float16 a, int n, float_status *status)
+{
+    FloatParts pa = float16_unpack_canonical(a, status);
+    FloatParts pr = scalbn_decomposed(pa, n, status);
+    return float16_round_pack_canonical(pr, status);
+}
+
+float32 float32_scalbn(float32 a, int n, float_status *status)
+{
+    FloatParts pa = float32_unpack_canonical(a, status);
+    FloatParts pr = scalbn_decomposed(pa, n, status);
+    return float32_round_pack_canonical(pr, status);
+}
+
+float64 float64_scalbn(float64 a, int n, float_status *status)
+{
+    FloatParts pa = float64_unpack_canonical(a, status);
+    FloatParts pr = scalbn_decomposed(pa, n, 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
@@ -6969,79 +7000,6 @@ MINMAX(32)
 MINMAX(64)
 
 
-/* Multiply A by 2 raised to the power N.  */
-float32 float32_scalbn(float32 a, int n, float_status *status)
-{
-    flag aSign;
-    int16_t aExp;
-    uint32_t aSig;
-
-    a = float32_squash_input_denormal(a, status);
-    aSig = extractFloat32Frac( a );
-    aExp = extractFloat32Exp( a );
-    aSign = extractFloat32Sign( a );
-
-    if ( aExp == 0xFF ) {
-        if ( aSig ) {
-            return propagateFloat32NaN(a, a, status);
-        }
-        return a;
-    }
-    if (aExp != 0) {
-        aSig |= 0x00800000;
-    } else if (aSig == 0) {
-        return a;
-    } else {
-        aExp++;
-    }
-
-    if (n > 0x200) {
-        n = 0x200;
-    } else if (n < -0x200) {
-        n = -0x200;
-    }
-
-    aExp += n - 1;
-    aSig <<= 7;
-    return normalizeRoundAndPackFloat32(aSign, aExp, aSig, status);
-}
-
-float64 float64_scalbn(float64 a, int n, float_status *status)
-{
-    flag aSign;
-    int16_t aExp;
-    uint64_t aSig;
-
-    a = float64_squash_input_denormal(a, status);
-    aSig = extractFloat64Frac( a );
-    aExp = extractFloat64Exp( a );
-    aSign = extractFloat64Sign( a );
-
-    if ( aExp == 0x7FF ) {
-        if ( aSig ) {
-            return propagateFloat64NaN(a, a, status);
-        }
-        return a;
-    }
-    if (aExp != 0) {
-        aSig |= LIT64( 0x0010000000000000 );
-    } else if (aSig == 0) {
-        return a;
-    } else {
-        aExp++;
-    }
-
-    if (n > 0x1000) {
-        n = 0x1000;
-    } else if (n < -0x1000) {
-        n = -0x1000;
-    }
-
-    aExp += n - 1;
-    aSig <<= 10;
-    return normalizeRoundAndPackFloat64(aSign, aExp, aSig, status);
-}
-
 floatx80 floatx80_scalbn(floatx80 a, int n, float_status *status)
 {
     flag aSign;
diff --git a/include/fpu/softfloat.h b/include/fpu/softfloat.h
index 3e6fdd756a..52621e0b79 100644
--- a/include/fpu/softfloat.h
+++ b/include/fpu/softfloat.h
@@ -244,6 +244,7 @@ float16 float16_sub(float16, float16, float_status *status);
 float16 float16_mul(float16, float16, float_status *status);
 float16 float16_muladd(float16, float16, float16, int, float_status *status);
 float16 float16_div(float16, float16, float_status *status);
+float16 float16_scalbn(float16, int, float_status *status);
 
 int float16_is_quiet_nan(float16, float_status *status);
 int float16_is_signaling_nan(float16, float_status *status);
-- 
2.15.1

On 01/24/2018 05:13 AM, Alex Bennée wrote:
> +/* Multiply A by 2 raised to the power N.  */
> +static FloatParts scalbn_decomposed(FloatParts a, int n,
> +                                          float_status *s)
> +{
> +    if (a.cls == float_class_normal) {
> +        a.exp += n;
> +    }
> +    return a;
> +}

...

> -    if ( aExp == 0x7FF ) {
> -        if ( aSig ) {
> -            return propagateFloat64NaN(a, a, status);

This is where we used to raise inexact.
In the new function we need to handle this as in round_to_int.


r~

Richard Henderson <richard.henderson@linaro.org> writes:

> On 01/24/2018 05:13 AM, Alex Bennée wrote:
>> +/* Multiply A by 2 raised to the power N.  */
>> +static FloatParts scalbn_decomposed(FloatParts a, int n,
>> +                                          float_status *s)
>> +{
>> +    if (a.cls == float_class_normal) {
>> +        a.exp += n;
>> +    }
>> +    return a;
>> +}
>
> ...
>
>> -    if ( aExp == 0x7FF ) {
>> -        if ( aSig ) {
>> -            return propagateFloat64NaN(a, a, status);
>
> This is where we used to raise inexact.
> In the new function we need to handle this as in round_to_int.

Why is the handling in round_canonical not good enough for this?

--
Alex Bennée

On 01/29/2018 03:59 AM, Alex Bennée wrote:
> 
> Richard Henderson <richard.henderson@linaro.org> writes:
> 
>> On 01/24/2018 05:13 AM, Alex Bennée wrote:
>>> +/* Multiply A by 2 raised to the power N.  */
>>> +static FloatParts scalbn_decomposed(FloatParts a, int n,
>>> +                                          float_status *s)
>>> +{
>>> +    if (a.cls == float_class_normal) {
>>> +        a.exp += n;
>>> +    }
>>> +    return a;
>>> +}
>>
>> ...
>>
>>> -    if ( aExp == 0x7FF ) {
>>> -        if ( aSig ) {
>>> -            return propagateFloat64NaN(a, a, status);
>>
>> This is where we used to raise inexact.
>> In the new function we need to handle this as in round_to_int.
> 
> Why is the handling in round_canonical not good enough for this?

This is scalbn -- you don't call round_canonical.


r~

On 01/29/2018 07:50 AM, Richard Henderson wrote:
> On 01/29/2018 03:59 AM, Alex Bennée wrote:
>>
>> Richard Henderson <richard.henderson@linaro.org> writes:
>>
>>> On 01/24/2018 05:13 AM, Alex Bennée wrote:
>>>> +/* Multiply A by 2 raised to the power N.  */
>>>> +static FloatParts scalbn_decomposed(FloatParts a, int n,
>>>> +                                          float_status *s)
>>>> +{
>>>> +    if (a.cls == float_class_normal) {
>>>> +        a.exp += n;
>>>> +    }
>>>> +    return a;
>>>> +}
>>>
>>> ...
>>>
>>>> -    if ( aExp == 0x7FF ) {
>>>> -        if ( aSig ) {
>>>> -            return propagateFloat64NaN(a, a, status);
>>>
>>> This is where we used to raise inexact.
>>> In the new function we need to handle this as in round_to_int.
>>
>> Why is the handling in round_canonical not good enough for this?
> 
> This is scalbn -- you don't call round_canonical.

Scratch that -- too early.  And I was trying to talk about invalid not inexact.
 What a cacophony of errors on my part.  No wonder we're not communicating.

But you'll notice that round_canonical does not raise invalid for nans.
Perhaps we could rearrange things so that it does, but at the moment it is the
responsibility of prior code to have set invalid and msnan/dnan.


r~