diff mbox

[v3,3/4] target-tilegx: Add double floating point implementation

Message ID 5669891F.6090707@emindsoft.com.cn
State New
Headers show

Commit Message

Chen Gang Dec. 10, 2015, 2:15 p.m. UTC
It passes gcc testsuite.

Signed-off-by: Chen Gang <gang.chen.5i5j@gmail.com>
---
 target-tilegx/helper-fdouble.c | 400 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 400 insertions(+)
 create mode 100644 target-tilegx/helper-fdouble.c

Comments

Richard Henderson Dec. 10, 2015, 9:17 p.m. UTC | #1
On 12/10/2015 06:15 AM, Chen Gang wrote:
> +#define TILEGX_F_MAN_HBIT   (1ULL << 59)
...
> +static uint64_t fr_to_man(float64 d)
> +{
> +    uint64_t val = get_f64_man(d) << 7;
> +
> +    if (get_f64_exp(d)) {
> +        val |= TILEGX_F_MAN_HBIT;
> +    }
> +
> +    return val;
> +}

One presumes that "HBIT" is the ieee implicit one bit.
A better name or better comments would help there.

Do we know for sure that "7" is the correct number of guard bits?  From the gcc
implementation of floatsidf, I might guess that the correct number is "4".

> +static uint32_t get_fdouble_vexp(uint64_t n)
> +{
> +    return extract32(n, 7, 13);
> +}

What's a "vexp"?

> +uint64_t helper_fdouble_unpack_min(CPUTLGState *env,
> +                                   uint64_t srca, uint64_t srcb)
> +{
> +    uint64_t v = 0;
> +    uint32_t expa = get_f64_exp(srca);
> +    uint32_t expb = get_f64_exp(srcb);
> +
> +    if (float64_is_any_nan(srca) || float64_is_any_nan(srcb)
> +        || float64_is_infinity(srca) || float64_is_infinity(srcb)) {
> +        return 0;
> +    } else if (expa > expb) {
> +        if (expa - expb < 64) {
> +            set_fdouble_man(&v, fr_to_man(srcb) >> (expa - expb));
> +        } else {
> +            return 0;
> +        }
> +    } else if (expa < expb) {
> +        if (expb - expa < 64) {
> +            set_fdouble_man(&v, fr_to_man(srca) >> (expb - expa));

I very sincerely doubt that a simple right-shift is correct.  In order to
obtain proper rounding for real computation, a sticky bit is required.  That
is, set bit 0 if any bits are shifted out.  See the implementation of
shift64RightJamming in fpu/softfloat-macros.h.

> +uint64_t helper_fdouble_addsub(CPUTLGState *env,
> +                               uint64_t dest, uint64_t srca, uint64_t srcb)
> +{
> +    if (get_fdouble_calc(srcb) == TILEGX_F_CALC_ADD) {
> +        return dest + srca; /* maybe set addsub overflow bit */

Definitely not.  That would be part of packing.

> +/* absolute-add/mul may cause add/mul carry or overflow */
> +static bool proc_oflow(uint64_t *flags, uint64_t *v, uint64_t *srcb)
> +{
> +    if (get_fdouble_man_of(*v)) {
> +        set_fdouble_vexp(flags, get_fdouble_vexp(*flags) + 1);
> +        *srcb >>= 1;
> +        *srcb |= *v << 63;
> +        *v >>= 1;
> +        clear_fdouble_man_of(v);
> +    }
> +    return get_fdouble_vexp(*flags) > TILEGX_F_EXP_DMAX;
> +}
> +
> +uint64_t helper_fdouble_pack2(CPUTLGState *env, uint64_t flags /* dest */,
> +                              uint64_t srca, uint64_t srcb)
> +{
> +    uint64_t v = srca;
> +    float64 d = float64_set_sign(float64_zero, get_fdouble_sign(flags));
> +
> +    /*
> +     * fdouble_add_flags, fdouble_sub_flags, or fdouble_mul_flags have
> +     * processed exceptions. So need not process fp_status, again.
> +     */

No need to process fp_status at all, actually.  Tile-GX (and pro) do not
support exception flags, so everything we do with fp_status is discarded.

Indeed, we should probably not store fp_status in env at all, but create it on
the stack in any function that actually needs one.


> +
> +    if (get_fdouble_nan(flags)) {
> +        return float64_val(float64_default_nan);
> +    } else if (get_fdouble_inf(flags)) {
> +        return float64_val(d |= float64_infinity);

s/|=/|/

> +    /* absolute-mul needs left shift 4 + 1 bytes to match the real mantissa */
> +    if (get_fdouble_calc(flags) == TILEGX_F_CALC_MUL) {
> +        v <<= 5;
> +        v |= srcb >> 59;
> +        srcb <<= 5;
> +    }

As with single, I don't like this calc thing.  We can infer what's required
from principals.

We're given two words containing mantissa, and a "flags" word containing sign,
exponent, and other flags.  For add, sub, and floatsidf, the compiler passes us
0 as the low word; for mul the compiler passes us the result of a 64x64->128
bit multiply.

The first step would be to normalize the 128-bit value so that the highest bit
set is TILEGX_F_MAN_HBIT in the high word, adjusting the exponent in the
process.  Fold the low word into the sticky bit of the high word (high |= (low
!= 0)) for rounding purposes.

The second step would be to round and pack, similar to roundAndPackFloat64,
except that your HBIT is at a different place than softfloat.c.

> +    d = calc(fsrca, fsrcb, fp_status); /* also check exceptions */

There are no exceptions to check.


r~
Chen Gang Dec. 11, 2015, 11:38 p.m. UTC | #2
On 12/11/15 05:17, Richard Henderson wrote:
> On 12/10/2015 06:15 AM, Chen Gang wrote:
>> +#define TILEGX_F_MAN_HBIT   (1ULL << 59)
> ...
>> +static uint64_t fr_to_man(float64 d)
>> +{
>> +    uint64_t val = get_f64_man(d) << 7;
>> +
>> +    if (get_f64_exp(d)) {
>> +        val |= TILEGX_F_MAN_HBIT;
>> +    }
>> +
>> +    return val;
>> +}
> 
> One presumes that "HBIT" is the ieee implicit one bit.
> A better name or better comments would help there.
> 

OK, thanks. And after think of again, I guess, the real hardware does
not use HBIT internally (use the full 64 bits as mantissa without HBIT).

But what I have done is still OK (use 59 bits + 1 HBIT as mantissa), for
59 bits are enough for double mantissa (52 bits). It makes the overflow
processing easier, but has to process mul operation specially.

It we have to try to match the real hardware have done, I shall rewrite
the related code for mantissa. (I guess, we need to match the real
hardware have done).


> Do we know for sure that "7" is the correct number of guard bits?  From the gcc
> implementation of floatsidf, I might guess that the correct number is "4".
> 

According to floatsidf, it seems "4", but after I expanded the bits, I
guess, it is "7".

/*
 * Double exp analyzing: (0x21b00 << 1) - 0x37(55) = 0x3ff
 *
 *   17  16  15  14  13  12  11  10   9   8   7    6   5   4   3   2   1   0
 *
 *    1   0   0   0   0   1   1   0   1   1   0    0   0   0   0   0   0   0
 *
 *    0   0   0   0   0   1   1   0   1   1   1    => 0x37(55)
 *
 *    0   1   1   1   1   1   1   1   1   1   1    => 0x3ff
 *
 */

I guess, I need restore this comments in helper_fdouble.c.


>> +static uint32_t get_fdouble_vexp(uint64_t n)
>> +{
>> +    return extract32(n, 7, 13);
>> +}
> 
> What's a "vexp"?
> 

It is exp + overflow bit + underflow bit. We can use vexp for internal
calculation, directly, and check uv and ov for the result. I guess the
real hardware will do like this.

The full description of the format is:

typedef union TileGXFPDFmtF {

    struct {
        uint64_t unknown0 : 7;    /* unknown */

        uint64_t vexp : 13;      /* vexp = exp | ov | uv */
#if 0 /* it is only the explanation for vexp above */
        uint64_t exp : 11;        /* exp, 0x21b << 1: 55 + TILEGX_F_EXP_DZERO */
        uint64_t ov : 1;          /* overflow for mul, low priority */
        uint64_t uv : 1;          /* underflow for mul, high priority */
#endif

        uint64_t sign : 1;        /* Sign bit for the total value */

        uint64_t calc: 2;         /* absolute add, sub, or mul */
        uint64_t inf: 1;          /* infinit */
        uint64_t nan: 1;          /* nan */

        /* Come from TILE-Gx ISA document, Table 7-2 for floating point */
        uint64_t unordered : 1;   /* The two are unordered */
        uint64_t lt : 1;          /* 1st is less than 2nd */
        uint64_t le : 1;          /* 1st is less than or equal to 2nd */
        uint64_t gt : 1;          /* 1st is greater than 2nd */
        uint64_t ge : 1;          /* 1st is greater than or equal to 2nd */
        uint64_t eq : 1;          /* The two operands are equal */
        uint64_t neq : 1;         /* The two operands are not equal */

        uint64_t unknown1 : 32;   /* unknown */
    } fmt;
    uint64_t ll;                  /* only for easy using */
} TileGXFPDFmtF;


>> +uint64_t helper_fdouble_unpack_min(CPUTLGState *env,
>> +                                   uint64_t srca, uint64_t srcb)
>> +{
>> +    uint64_t v = 0;
>> +    uint32_t expa = get_f64_exp(srca);
>> +    uint32_t expb = get_f64_exp(srcb);
>> +
>> +    if (float64_is_any_nan(srca) || float64_is_any_nan(srcb)
>> +        || float64_is_infinity(srca) || float64_is_infinity(srcb)) {
>> +        return 0;
>> +    } else if (expa > expb) {
>> +        if (expa - expb < 64) {
>> +            set_fdouble_man(&v, fr_to_man(srcb) >> (expa - expb));
>> +        } else {
>> +            return 0;
>> +        }
>> +    } else if (expa < expb) {
>> +        if (expb - expa < 64) {
>> +            set_fdouble_man(&v, fr_to_man(srca) >> (expb - expa));
> 
> I very sincerely doubt that a simple right-shift is correct.  In order to
> obtain proper rounding for real computation, a sticky bit is required.  That
> is, set bit 0 if any bits are shifted out.  See the implementation of
> shift64RightJamming in fpu/softfloat-macros.h.
> 

Oh, really, thanks.


>> +uint64_t helper_fdouble_addsub(CPUTLGState *env,
>> +                               uint64_t dest, uint64_t srca, uint64_t srcb)
>> +{
>> +    if (get_fdouble_calc(srcb) == TILEGX_F_CALC_ADD) {
>> +        return dest + srca; /* maybe set addsub overflow bit */
> 
> Definitely not.  That would be part of packing.
> 

If we need to try to match the real hardware have done, the related
implementation above is incorrect.

And for my current implementation (I guess, it should be correct):

typedef union TileGXFPDFmtV {
    struct {
        uint64_t mantissa : 60;   /* mantissa */
        uint64_t overflow : 1;    /* carry/overflow bit for absolute add/mul */
        uint64_t unknown1 : 3;    /* unknown */
    } fmt;
    uint64_t ll;                  /* only for easy using */
} TileGXFPDFmtV;


In helper_fdouble_addsub(), both dest and srca are unpacked, so they are
within 60 bits. So one time absolute add are within 61 bits, so let bit
61 as overflow bit is enough.


>> +/* absolute-add/mul may cause add/mul carry or overflow */
>> +static bool proc_oflow(uint64_t *flags, uint64_t *v, uint64_t *srcb)
>> +{
>> +    if (get_fdouble_man_of(*v)) {
>> +        set_fdouble_vexp(flags, get_fdouble_vexp(*flags) + 1);
>> +        *srcb >>= 1;
>> +        *srcb |= *v << 63;
>> +        *v >>= 1;
>> +        clear_fdouble_man_of(v);
>> +    }
>> +    return get_fdouble_vexp(*flags) > TILEGX_F_EXP_DMAX;
>> +}
>> +
>> +uint64_t helper_fdouble_pack2(CPUTLGState *env, uint64_t flags /* dest */,
>> +                              uint64_t srca, uint64_t srcb)
>> +{
>> +    uint64_t v = srca;
>> +    float64 d = float64_set_sign(float64_zero, get_fdouble_sign(flags));
>> +
>> +    /*
>> +     * fdouble_add_flags, fdouble_sub_flags, or fdouble_mul_flags have
>> +     * processed exceptions. So need not process fp_status, again.
>> +     */
> 
> No need to process fp_status at all, actually.  Tile-GX (and pro) do not
> support exception flags, so everything we do with fp_status is discarded.
> 
> Indeed, we should probably not store fp_status in env at all, but create it on
> the stack in any function that actually needs one.
> 

OK, thanks.

>> +
>> +    if (get_fdouble_nan(flags)) {
>> +        return float64_val(float64_default_nan);
>> +    } else if (get_fdouble_inf(flags)) {
>> +        return float64_val(d |= float64_infinity);
> 
> s/|=/|/
> 

OK, thanks.

>> +    /* absolute-mul needs left shift 4 + 1 bytes to match the real mantissa */
>> +    if (get_fdouble_calc(flags) == TILEGX_F_CALC_MUL) {
>> +        v <<= 5;
>> +        v |= srcb >> 59;
>> +        srcb <<= 5;
>> +    }
> 
> As with single, I don't like this calc thing.  We can infer what's required
> from principals.
> 
> We're given two words containing mantissa, and a "flags" word containing sign,
> exponent, and other flags.  For add, sub, and floatsidf, the compiler passes us
> 0 as the low word; for mul the compiler passes us the result of a 64x64->128
> bit multiply.
> 

OK, thanks. It looks, we have to try to match what the hardware have
done.

> The first step would be to normalize the 128-bit value so that the highest bit
> set is TILEGX_F_MAN_HBIT in the high word, adjusting the exponent in the
> process.  Fold the low word into the sticky bit of the high word (high |= (low
> != 0)) for rounding purposes.
> 

OK, thanks. And my original implementation did not consider about the
sticky bit.

> The second step would be to round and pack, similar to roundAndPackFloat64,
> except that your HBIT is at a different place than softfloat.c.
> 

It sounds good (and originally I really considered about it). If we have
an export common function for it, that will be really good.

At present, I use (u)int64_to_float64(), then process exp again.


>> +    d = calc(fsrca, fsrcb, fp_status); /* also check exceptions */
> 
> There are no exceptions to check.
>

OK, thanks.
Richard Henderson Dec. 12, 2015, 12:41 a.m. UTC | #3
On 12/11/2015 03:38 PM, Chen Gang wrote:
>
> On 12/11/15 05:17, Richard Henderson wrote:
>> On 12/10/2015 06:15 AM, Chen Gang wrote:
>>> +#define TILEGX_F_MAN_HBIT   (1ULL << 59)
>> ...
>>> +static uint64_t fr_to_man(float64 d)
>>> +{
>>> +    uint64_t val = get_f64_man(d) << 7;
>>> +
>>> +    if (get_f64_exp(d)) {
>>> +        val |= TILEGX_F_MAN_HBIT;
>>> +    }
>>> +
>>> +    return val;
>>> +}
>>
>> One presumes that "HBIT" is the ieee implicit one bit.
>> A better name or better comments would help there.
>>
>
> OK, thanks. And after think of again, I guess, the real hardware does
> not use HBIT internally (use the full 64 bits as mantissa without HBIT).

It must do.  Otherwise the arithmetic doesn't work out.

> But what I have done is still OK (use 59 bits + 1 HBIT as mantissa), for
> 59 bits are enough for double mantissa (52 bits). It makes the overflow
> processing easier, but has to process mul operation specially.

What you have works.  But the mul operation isn't as special as you make it out 
-- aside from requiring at least 104 bits as intermediate -- in that when one 
implements what the hardware does, subtraction also may require significant 
normalization.

> According to floatsidf, it seems "4", but after I expanded the bits, I
> guess, it is "7".
>
> /*
>   * Double exp analyzing: (0x21b00 << 1) - 0x37(55) = 0x3ff
>   *
>   *   17  16  15  14  13  12  11  10   9   8   7    6   5   4   3   2   1   0
>   *
>   *    1   0   0   0   0   1   1   0   1   1   0    0   0   0   0   0   0   0
>   *
>   *    0   0   0   0   0   1   1   0   1   1   1    => 0x37(55)
>   *
>   *    0   1   1   1   1   1   1   1   1   1   1    => 0x3ff
>   *
>   */

That's the exponent within the flags temporary.  It has nothing to do with the 
position of the extracted mantissa.

FWIW, the minimum shift would be 3, in order to properly implement rounding; if 
the hardware uses a shift of 4, that's fine too.

What I would love to know is if the shift present in floatsidf is not really 
required; equally valid to adjust 0x21b00 by 4.  Meaning normalization would do 
a proper job with the entire given mantissa.  This would require better 
documentation, or access to hardware to verify.

>>> +uint64_t helper_fdouble_addsub(CPUTLGState
> And for my current implementation (I guess, it should be correct):
>
> typedef union TileGXFPDFmtV {
>      struct {
>          uint64_t mantissa : 60;   /* mantissa */
>          uint64_t overflow : 1;    /* carry/overflow bit for absolute add/mul */
>          uint64_t unknown1 : 3;    /* unknown */

I personally like to call all 4 of the top bits overflow.  But I have no idea 
what the real hardware actually does.

> In helper_fdouble_addsub(), both dest and srca are unpacked, so they are
> within 60 bits. So one time absolute add are within 61 bits, so let bit
> 61 as overflow bit is enough.

True.  But if all 4 top bits are considered overflow, then one could implement 
floatdidf fairly easily.  But I suspect that real hw doesn't work that way, or 
it would have already been done.


r~
Chen Gang Dec. 12, 2015, 2:45 a.m. UTC | #4
On 12/12/15 08:41, Richard Henderson wrote:
> On 12/11/2015 03:38 PM, Chen Gang wrote:
>>
>> On 12/11/15 05:17, Richard Henderson wrote:
>>> On 12/10/2015 06:15 AM, Chen Gang wrote:
>>>> +#define TILEGX_F_MAN_HBIT   (1ULL << 59)
>>> ...
>>>> +static uint64_t fr_to_man(float64 d)
>>>> +{
>>>> +    uint64_t val = get_f64_man(d) << 7;
>>>> +
>>>> +    if (get_f64_exp(d)) {
>>>> +        val |= TILEGX_F_MAN_HBIT;
>>>> +    }
>>>> +
>>>> +    return val;
>>>> +}
>>>
>>> One presumes that "HBIT" is the ieee implicit one bit.
>>> A better name or better comments would help there.
>>>
>>
>> OK, thanks. And after think of again, I guess, the real hardware does
>> not use HBIT internally (use the full 64 bits as mantissa without HBIT).
> 
> It must do.  Otherwise the arithmetic doesn't work out.
> 

Oh, yes, and we have to use my original implementation (60 for mantissa,
4 bits for other using).

>> But what I have done is still OK (use 59 bits + 1 HBIT as mantissa), for
>> 59 bits are enough for double mantissa (52 bits). It makes the overflow
>> processing easier, but has to process mul operation specially.
> 
> What you have works.  But the mul operation isn't as special as you make it out -- aside from requiring at least 104 bits as intermediate -- in that when one implements what the hardware does, subtraction also may require significant normalization.
> 

I guess, you misunderstood what I said (my English is not quite well).

For mul, at least, it needs (104 - 1) bits, At present, we have 120 bits
for it (in fact, our mul generates 119 bits result). So it is enough.


>> According to floatsidf, it seems "4", but after I expanded the bits, I
>> guess, it is "7".
>>
>> /*
>>   * Double exp analyzing: (0x21b00 << 1) - 0x37(55) = 0x3ff
>>   *
>>   *   17  16  15  14  13  12  11  10   9   8   7    6   5   4   3   2   1   0
>>   *
>>   *    1   0   0   0   0   1   1   0   1   1   0    0   0   0   0   0   0   0
>>   *
>>   *    0   0   0   0   0   1   1   0   1   1   1    => 0x37(55)
>>   *
>>   *    0   1   1   1   1   1   1   1   1   1   1    => 0x3ff
>>   *
>>   */
> 
> That's the exponent within the flags temporary.  It has nothing to do with the position of the extracted mantissa.
> 

0x37(55) + 4 (guard bits) + 1 (HBIT) = 60 bits.

So, if the above is correct, the mantissa is 60 bits (with HBIT), and
bit 18 in flags for overflow, bit 19 for underflow (bit 20 must be for
sign).

> FWIW, the minimum shift would be 3, in order to properly implement rounding; if the hardware uses a shift of 4, that's fine too.
> 

I guess, so it uses 4 guard bits.

> What I would love to know is if the shift present in floatsidf is not really required; equally valid to adjust 0x21b00 by 4.  Meaning normalization would do a proper job with the entire given mantissa.  This would require better documentation, or access to hardware to verify.
> 

I guess, before call any fdouble insns, we can use the low 4 bits as
mantissa (e.g. calc mul), but when call any fdouble insn, we can not use
the lower 4 guard bits, so floatsidf has to shift 4 bits left.

>>>> +uint64_t helper_fdouble_addsub(CPUTLGState
>> And for my current implementation (I guess, it should be correct):
>>
>> typedef union TileGXFPDFmtV {
>>      struct {
>>          uint64_t mantissa : 60;   /* mantissa */
>>          uint64_t overflow : 1;    /* carry/overflow bit for absolute add/mul */
>>          uint64_t unknown1 : 3;    /* unknown */
> 
> I personally like to call all 4 of the top bits overflow.  But I have no idea what the real hardware actually does.
> 
>> In helper_fdouble_addsub(), both dest and srca are unpacked, so they are
>> within 60 bits. So one time absolute add are within 61 bits, so let bit
>> 61 as overflow bit is enough.
> 
> True.  But if all 4 top bits are considered overflow, then one could implement floatdidf fairly easily.  But I suspect that real hw doesn't work that way, or it would have already been done.
> 

So, I only assumed bit 60 is for overflow, the high 3 bits are unknown.

For me, if one bit for overflow is enough, the hardware will save the
other bits for another using (or are reserved for future).


Thanks.
diff mbox

Patch

diff --git a/target-tilegx/helper-fdouble.c b/target-tilegx/helper-fdouble.c
new file mode 100644
index 0000000..3b824f7
--- /dev/null
+++ b/target-tilegx/helper-fdouble.c
@@ -0,0 +1,400 @@ 
+/*
+ * QEMU TILE-Gx helpers
+ *
+ *  Copyright (c) 2015 Chen Gang
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, see
+ * <http://www.gnu.org/licenses/lgpl-2.1.html>
+ */
+
+#include "cpu.h"
+#include "qemu-common.h"
+#include "exec/helper-proto.h"
+#include "fpu/softfloat.h"
+
+#include "helper-fshared.c"
+
+/*
+ * FDouble instructions implemenation:
+ *
+ * fdouble_unpack_min   ; srca and srcb are float_64 value.
+ *                      ; get the min absolute value's mantissa.
+ *                      ; move "mantissa >> (exp_max - exp_min)" to dest.
+ *
+ * fdouble_unpack_max   ; srca and srcb are float_64 value.
+ *                      ; get the max absolute value's mantissa.
+ *                      ; move mantissa to dest.
+ *
+ * fdouble_add_flags    ; srca and srcb are float_64 value.
+ *                      ; calc exp (exp_max), sign, and comp bits for flags.
+ *                      ; set addsub bit to flags and move flags to dest.
+ *
+ * fdouble_sub_flags    ; srca and srcb are float_64 value.
+ *                      ; calc exp (exp_max), sign, and comp bits for flags.
+ *                      ; set addsub bit to flags and move flags to dest.
+ *
+ * fdouble_addsub:      ; dest, srca (max, min mantissa), and srcb (flags).
+ *                      ; "dest +/- srca" depend on the add/sub bit of flags.
+ *                      ; move result mantissa to dest.
+ *
+ * fdouble_mul_flags:   ; srca and srcb are float_64 value.
+ *                      ; calc sign (xor), exp (min + max), and comp bits.
+ *                      ; mix sign, exp, and comp bits as flags to dest.
+ *
+ * fdouble_pack1        ; move srcb (flags) to dest.
+ *
+ * fdouble_pack2        ; srca, srcb (high, low mantissa), and dest (flags)
+ *                      ; normalize and pack result from srca, srcb, and dest.
+ *                      ; move result to dest.
+ */
+
+#define TILEGX_F_EXP_DZERO  0x3ff /* Zero exp for double 11-bits */
+#define TILEGX_F_EXP_DMAX   0x7fe /* max exp for double 11-bits */
+#define TILEGX_F_EXP_DUF    0x1000/* underflow exp bit for double */
+
+#define TILEGX_F_MAN_HBIT   (1ULL << 59)
+
+#define TILEGX_F_CALC_ADD   1     /* Perform absolute add operation */
+#define TILEGX_F_CALC_SUB   2     /* Perform absolute sub operation */
+#define TILEGX_F_CALC_MUL   3     /* Perform absolute mul operation */
+
+static uint32_t get_f64_exp(float64 d)
+{
+    return extract64(float64_val(d), 52, 11);
+}
+
+static void set_f64_exp(float64 *d, uint32_t exp)
+{
+    *d = make_float64(deposit64(float64_val(*d), 52, 11, exp));
+}
+
+static uint64_t get_f64_man(float64 d)
+{
+    return extract64(float64_val(d), 0, 52);
+}
+
+static uint64_t fr_to_man(float64 d)
+{
+    uint64_t val = get_f64_man(d) << 7;
+
+    if (get_f64_exp(d)) {
+        val |= TILEGX_F_MAN_HBIT;
+    }
+
+    return val;
+}
+
+static uint64_t get_fdouble_man(uint64_t n)
+{
+    return extract64(n, 0, 60);
+}
+
+static void set_fdouble_man(uint64_t *n, uint64_t man)
+{
+    *n = deposit64(*n, 0, 60, man);
+}
+
+static uint64_t get_fdouble_man_of(uint64_t n)
+{
+    return test_bit(60, &n);
+}
+
+static void clear_fdouble_man_of(uint64_t *n)
+{
+    return clear_bit(60, n);
+}
+
+static uint32_t get_fdouble_nan(uint64_t n)
+{
+    return test_bit(24, &n);
+}
+
+static void set_fdouble_nan(uint64_t *n)
+{
+    set_bit(24, n);
+}
+
+static uint32_t get_fdouble_inf(uint64_t n)
+{
+    return test_bit(23, &n);
+}
+
+static void set_fdouble_inf(uint64_t *n)
+{
+    set_bit(23, n);
+}
+
+static uint32_t get_fdouble_calc(uint64_t n)
+{
+    return extract32(n, 21, 2);
+}
+
+static void set_fdouble_calc(uint64_t *n, uint32_t calc)
+{
+    *n = deposit64(*n, 21, 2, calc);
+}
+
+static uint32_t get_fdouble_sign(uint64_t n)
+{
+    return test_bit(20, &n);
+}
+
+static void set_fdouble_sign(uint64_t *n)
+{
+    set_bit(20, n);
+}
+
+static uint32_t get_fdouble_vexp(uint64_t n)
+{
+    return extract32(n, 7, 13);
+}
+
+static void set_fdouble_vexp(uint64_t *n, uint32_t vexp)
+{
+    *n = deposit64(*n, 7, 13, vexp);
+}
+
+uint64_t helper_fdouble_unpack_min(CPUTLGState *env,
+                                   uint64_t srca, uint64_t srcb)
+{
+    uint64_t v = 0;
+    uint32_t expa = get_f64_exp(srca);
+    uint32_t expb = get_f64_exp(srcb);
+
+    if (float64_is_any_nan(srca) || float64_is_any_nan(srcb)
+        || float64_is_infinity(srca) || float64_is_infinity(srcb)) {
+        return 0;
+    } else if (expa > expb) {
+        if (expa - expb < 64) {
+            set_fdouble_man(&v, fr_to_man(srcb) >> (expa - expb));
+        } else {
+            return 0;
+        }
+    } else if (expa < expb) {
+        if (expb - expa < 64) {
+            set_fdouble_man(&v, fr_to_man(srca) >> (expb - expa));
+        } else {
+            return 0;
+        }
+    } else if (get_f64_man(srca) > get_f64_man(srcb)) {
+        set_fdouble_man(&v, fr_to_man(srcb));
+    } else {
+        set_fdouble_man(&v, fr_to_man(srca));
+    }
+
+    return v;
+}
+
+uint64_t helper_fdouble_unpack_max(CPUTLGState *env,
+                                   uint64_t srca, uint64_t srcb)
+{
+    uint64_t v = 0;
+    uint32_t expa = get_f64_exp(srca);
+    uint32_t expb = get_f64_exp(srcb);
+
+    if (float64_is_any_nan(srca) || float64_is_any_nan(srcb)
+        || float64_is_infinity(srca) || float64_is_infinity(srcb)) {
+        return 0;
+    } else if (expa > expb) {
+        set_fdouble_man(&v, fr_to_man(srca));
+    } else if (expa < expb) {
+        set_fdouble_man(&v, fr_to_man(srcb));
+    } else if (get_f64_man(srca) > get_f64_man(srcb)) {
+        set_fdouble_man(&v, fr_to_man(srca));
+    } else {
+        set_fdouble_man(&v, fr_to_man(srcb));
+    }
+
+    return v;
+}
+
+uint64_t helper_fdouble_addsub(CPUTLGState *env,
+                               uint64_t dest, uint64_t srca, uint64_t srcb)
+{
+    if (get_fdouble_calc(srcb) == TILEGX_F_CALC_ADD) {
+        return dest + srca; /* maybe set addsub overflow bit */
+    } else {
+        return dest - srca;
+    }
+}
+
+/* absolute-add/mul may cause add/mul carry or overflow */
+static bool proc_oflow(uint64_t *flags, uint64_t *v, uint64_t *srcb)
+{
+    if (get_fdouble_man_of(*v)) {
+        set_fdouble_vexp(flags, get_fdouble_vexp(*flags) + 1);
+        *srcb >>= 1;
+        *srcb |= *v << 63;
+        *v >>= 1;
+        clear_fdouble_man_of(v);
+    }
+    return get_fdouble_vexp(*flags) > TILEGX_F_EXP_DMAX;
+}
+
+uint64_t helper_fdouble_pack2(CPUTLGState *env, uint64_t flags /* dest */,
+                              uint64_t srca, uint64_t srcb)
+{
+    uint64_t v = srca;
+    float64 d = float64_set_sign(float64_zero, get_fdouble_sign(flags));
+
+    /*
+     * fdouble_add_flags, fdouble_sub_flags, or fdouble_mul_flags have
+     * processed exceptions. So need not process fp_status, again.
+     */
+
+    if (get_fdouble_nan(flags)) {
+        return float64_val(float64_default_nan);
+    } else if (get_fdouble_inf(flags)) {
+        return float64_val(d |= float64_infinity);
+    }
+
+    /* absolute-mul needs left shift 4 + 1 bytes to match the real mantissa */
+    if (get_fdouble_calc(flags) == TILEGX_F_CALC_MUL) {
+        v <<= 5;
+        v |= srcb >> 59;
+        srcb <<= 5;
+    }
+
+    /* must check underflow, firstly */
+    if (get_fdouble_vexp(flags) & TILEGX_F_EXP_DUF) {
+        return float64_val(d);
+    }
+
+    if (proc_oflow(&flags, &v, &srcb)) {
+        return float64_val(d |= float64_infinity);
+    }
+
+    while (!(get_fdouble_man(v) & TILEGX_F_MAN_HBIT)
+           && (get_fdouble_man(v) | srcb)) {
+        set_fdouble_vexp(&flags, get_fdouble_vexp(flags) - 1);
+        set_fdouble_man(&v, get_fdouble_man(v) << 1);
+        set_fdouble_man(&v, get_fdouble_man(v) | (srcb >> 63));
+        srcb <<= 1;
+    }
+
+    /* check underflow, again, after format */
+    if ((get_fdouble_vexp(flags) & TILEGX_F_EXP_DUF) || !get_fdouble_man(v)) {
+        return float64_val(d);
+    }
+
+    if (get_fdouble_sign(flags)) {
+        d = int64_to_float64(0 - get_fdouble_man(v), &env->fp_status);
+    } else {
+        d = uint64_to_float64(get_fdouble_man(v), &env->fp_status);
+    }
+
+    if (get_f64_exp(d) == 59 + TILEGX_F_EXP_DZERO) {
+        set_f64_exp(&d, get_fdouble_vexp(flags));
+    } else {                            /* for carry and overflow again */
+        set_f64_exp(&d, get_fdouble_vexp(flags) + 1);
+        if (get_f64_exp(d) == TILEGX_F_EXP_DMAX) {
+            d = float64_infinity;
+        }
+    }
+
+    d = float64_set_sign(d, get_fdouble_sign(flags));
+
+    return float64_val(d);
+}
+
+static void ana_bits(float_status *fp_status,
+                     float64 fsrca, float64 fsrcb, uint64_t *dfmt)
+{
+    if (float64_eq(fsrca, fsrcb, fp_status)) {
+        *dfmt |= create_fsfd_flag_eq();
+    } else {
+        *dfmt |= create_fsfd_flag_ne();
+    }
+
+    if (float64_lt(fsrca, fsrcb, fp_status)) {
+        *dfmt |= create_fsfd_flag_lt();
+    }
+    if (float64_le(fsrca, fsrcb, fp_status)) {
+        *dfmt |= create_fsfd_flag_le();
+    }
+
+    if (float64_lt(fsrcb, fsrca, fp_status)) {
+        *dfmt |= create_fsfd_flag_gt();
+    }
+    if (float64_le(fsrcb, fsrca, fp_status)) {
+        *dfmt |= create_fsfd_flag_ge();
+    }
+
+    if (float64_unordered(fsrca, fsrcb, fp_status)) {
+        *dfmt |= create_fsfd_flag_un();
+    }
+}
+
+static uint64_t main_calc(float_status *fp_status,
+                          float64 fsrca, float64 fsrcb,
+                          float64 (*calc)(float64, float64, float_status *))
+{
+    float64 d;
+    uint64_t flags = 0;
+    uint32_t expa = get_f64_exp(fsrca);
+    uint32_t expb = get_f64_exp(fsrcb);
+
+    ana_bits(fp_status, fsrca, fsrcb, &flags);
+
+    d = calc(fsrca, fsrcb, fp_status); /* also check exceptions */
+    if (float64_is_neg(d)) {
+        set_fdouble_sign(&flags);
+    }
+
+    if (float64_is_any_nan(d)) {
+        set_fdouble_nan(&flags);
+    } else if (float64_is_infinity(d)) {
+        set_fdouble_inf(&flags);
+    } else if (calc == float64_add) {
+        set_fdouble_vexp(&flags, (expa > expb) ? expa : expb);
+        set_fdouble_calc(&flags,
+                         (float64_is_neg(fsrca) == float64_is_neg(fsrcb))
+                             ? TILEGX_F_CALC_ADD : TILEGX_F_CALC_SUB);
+
+    } else if (calc == float64_sub) {
+        set_fdouble_vexp(&flags, (expa > expb) ? expa : expb);
+        set_fdouble_calc(&flags,
+                         (float64_is_neg(fsrca) != float64_is_neg(fsrcb))
+                             ? TILEGX_F_CALC_ADD : TILEGX_F_CALC_SUB);
+
+    } else {
+        set_fdouble_vexp(&flags, (int64_t)(expa - TILEGX_F_EXP_DZERO)
+                                 + (int64_t)(expb - TILEGX_F_EXP_DZERO)
+                                 + TILEGX_F_EXP_DZERO);
+        set_fdouble_calc(&flags, TILEGX_F_CALC_MUL);
+    }
+
+    return flags;
+}
+
+uint64_t helper_fdouble_add_flags(CPUTLGState *env,
+                                  uint64_t srca, uint64_t srcb)
+{
+    return main_calc(&env->fp_status,
+                     make_float64(srca), make_float64(srcb), float64_add);
+}
+
+uint64_t helper_fdouble_sub_flags(CPUTLGState *env,
+                                  uint64_t srca, uint64_t srcb)
+{
+    return main_calc(&env->fp_status,
+                     make_float64(srca), make_float64(srcb), float64_sub);
+}
+
+uint64_t helper_fdouble_mul_flags(CPUTLGState *env,
+                                  uint64_t srca, uint64_t srcb)
+{
+    return main_calc(&env->fp_status,
+                     make_float64(srca), make_float64(srcb), float64_mul);
+}