Re: [Qemu-devel] [PATCH] powerpc: fix denorm float->double conversion

Re: [Qemu-devel] [PATCH] powerpc: fix denorm float->double conversion

[David Gibson]

[-- Attachment #1: Type: text/plain, Size: 3901 bytes --]

On Sat, Mar 23, 2019 at 10:24:11PM +0000, Sergei Trofimovich wrote:
> The bug is initially discovered in GHC test suite. Here is minimal reproducer:
> 
> ```c
> 
> int main() {
>     volatile float f;
>     volatile double d;
> 
>     *(volatile uint32_t*)&f = 0xc0de;
>     d = f;
>     printf("f  = %#x\n", *(volatile uint32_t*)&f);
>     printf("d  = %#llx (expect 0x37981bc000000000)\n",
>         *(volatile uint64_t*)&d);
>     printf("d  = %e\n", d);
>     f = d;
>     printf("f  = %#x\n", *(volatile uint32_t*)&f);
> }
> ```
> 
> ```
> $ powerpc-unknown-linux-gnu-gcc -O2 a.c -Wall -o a \
>     -fno-strict-aliasing -static && qemu-ppc ./a
> f  = 0xc0de
> d  = 0x37a00000000c0de0 (expect 0x37981bc000000000)
> d  = 9.183550e-41
> f  = 0x10000
> ```
> 
> Here denormalization conversion has a few bugs:
> - significand (abs_arg) has 32-bit unsigned wraparound in
>     ret |= abs_arg << (shift + 29);
> - significand does not drop explicit leading '1' in denorm
>   'float' when converting to normalized 'double'
> - significand had an off-by-one shift
> 
> CC: Richard Henderson <richard.henderson@linaro.org>
> CC: David Gibson <david@gibson.dropbear.id.au>
> CC: qemu-ppc@nongnu.org
> CC: qemu-devel@nongnu.org
> Bug: https://bugs.launchpad.net/qemu/+bug/1821444
> Signed-off-by: Sergei Trofimovich <slyfox@gentoo.org>

LGTM, but I don't know much about floating point.

Richard, can you review this?

> ---
>  target/ppc/fpu_helper.c | 32 +++++++++++++++++++++++++++-----
>  1 file changed, 27 insertions(+), 5 deletions(-)
> 
> diff --git a/target/ppc/fpu_helper.c b/target/ppc/fpu_helper.c
> index 2ed4f42275..1e8b014890 100644
> --- a/target/ppc/fpu_helper.c
> +++ b/target/ppc/fpu_helper.c
> @@ -64,13 +64,35 @@ uint64_t helper_todouble(uint32_t arg)
>          ret |= (uint64_t)extract32(arg, 0, 30) << 29;
>      } else {
>          /* Zero or Denormalized operand.  */
> -        ret = (uint64_t)extract32(arg, 31, 1) << 63;
> +
> +        /*
> +         * Conversion mechanics:
> +         * float denorm (2^(-126) - biased):
> +         *    [ sign (1 bit) | exp32 (8 bits)  | sign32 (23 bits) ]
> +         *                 s                0    0001abc...def
> +         * double norm (2^(-1023) - biased):
> +         *    [ sign (1 bit) | exp64 (11 bits) | sign64 (52 bits) ]
> +         *                 s              exp    abc...def 00..0
> +         * Thus we are performing the following conversion steps:
> +         * 1. preserve the sign
> +         * 2. normalize denorm sign32:
> +         *   2a. drop explicit leading '1' as normalized numbers
> +         *       don't contain it
> +         *   2b. calculate the bit-shift needed to match implicit '1'
> +         * 3. calculate 'exp64' as bias delta plus denorm offset
> +         * 4. put calculated 'sign64' into new location
> +         */
> +        ret = (uint64_t)extract32(arg, 31, 1) << 63; /* [1.] */
>          if (unlikely(abs_arg != 0)) {
>              /* Denormalized operand.  */
> -            int shift = clz32(abs_arg) - 9;
> -            int exp = -126 - shift + 1023;
> -            ret |= (uint64_t)exp << 52;
> -            ret |= abs_arg << (shift + 29);
> +            int lz = clz32(abs_arg);
> +            abs_arg &= ~(1 << (31 - lz)); /* [2a.] */
> +
> +            /* shift within sign32 includeing leading '1' */
> +            int shift = lz + 1 - (32 - 23);
> +            int exp = -126 + 1023 - shift; /* [2b]. */
> +            ret |= (uint64_t)exp << 52; /* [3.] */
> +            ret |= (uint64_t)abs_arg << (52 - 23 + shift); /* [4.] */
>          }
>      }
>      return ret;

-- 
David Gibson			| I'll have my music baroque, and my code
david AT gibson.dropbear.id.au	| minimalist, thank you.  NOT _the_ _other_
				| _way_ _around_!
http://www.ozlabs.org/~dgibson

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 833 bytes --]

Re: [Qemu-devel] [PATCH] powerpc: fix denorm float->double conversion

[David Gibson]

[-- Attachment #1: Type: text/plain, Size: 3901 bytes --]

On Sat, Mar 23, 2019 at 10:24:11PM +0000, Sergei Trofimovich wrote:
> The bug is initially discovered in GHC test suite. Here is minimal reproducer:
> 
> ```c
> 
> int main() {
>     volatile float f;
>     volatile double d;
> 
>     *(volatile uint32_t*)&f = 0xc0de;
>     d = f;
>     printf("f  = %#x\n", *(volatile uint32_t*)&f);
>     printf("d  = %#llx (expect 0x37981bc000000000)\n",
>         *(volatile uint64_t*)&d);
>     printf("d  = %e\n", d);
>     f = d;
>     printf("f  = %#x\n", *(volatile uint32_t*)&f);
> }
> ```
> 
> ```
> $ powerpc-unknown-linux-gnu-gcc -O2 a.c -Wall -o a \
>     -fno-strict-aliasing -static && qemu-ppc ./a
> f  = 0xc0de
> d  = 0x37a00000000c0de0 (expect 0x37981bc000000000)
> d  = 9.183550e-41
> f  = 0x10000
> ```
> 
> Here denormalization conversion has a few bugs:
> - significand (abs_arg) has 32-bit unsigned wraparound in
>     ret |= abs_arg << (shift + 29);
> - significand does not drop explicit leading '1' in denorm
>   'float' when converting to normalized 'double'
> - significand had an off-by-one shift
> 
> CC: Richard Henderson <richard.henderson@linaro.org>
> CC: David Gibson <david@gibson.dropbear.id.au>
> CC: qemu-ppc@nongnu.org
> CC: qemu-devel@nongnu.org
> Bug: https://bugs.launchpad.net/qemu/+bug/1821444
> Signed-off-by: Sergei Trofimovich <slyfox@gentoo.org>

LGTM, but I don't know much about floating point.

Richard, can you review this?

> ---
>  target/ppc/fpu_helper.c | 32 +++++++++++++++++++++++++++-----
>  1 file changed, 27 insertions(+), 5 deletions(-)
> 
> diff --git a/target/ppc/fpu_helper.c b/target/ppc/fpu_helper.c
> index 2ed4f42275..1e8b014890 100644
> --- a/target/ppc/fpu_helper.c
> +++ b/target/ppc/fpu_helper.c
> @@ -64,13 +64,35 @@ uint64_t helper_todouble(uint32_t arg)
>          ret |= (uint64_t)extract32(arg, 0, 30) << 29;
>      } else {
>          /* Zero or Denormalized operand.  */
> -        ret = (uint64_t)extract32(arg, 31, 1) << 63;
> +
> +        /*
> +         * Conversion mechanics:
> +         * float denorm (2^(-126) - biased):
> +         *    [ sign (1 bit) | exp32 (8 bits)  | sign32 (23 bits) ]
> +         *                 s                0    0001abc...def
> +         * double norm (2^(-1023) - biased):
> +         *    [ sign (1 bit) | exp64 (11 bits) | sign64 (52 bits) ]
> +         *                 s              exp    abc...def 00..0
> +         * Thus we are performing the following conversion steps:
> +         * 1. preserve the sign
> +         * 2. normalize denorm sign32:
> +         *   2a. drop explicit leading '1' as normalized numbers
> +         *       don't contain it
> +         *   2b. calculate the bit-shift needed to match implicit '1'
> +         * 3. calculate 'exp64' as bias delta plus denorm offset
> +         * 4. put calculated 'sign64' into new location
> +         */
> +        ret = (uint64_t)extract32(arg, 31, 1) << 63; /* [1.] */
>          if (unlikely(abs_arg != 0)) {
>              /* Denormalized operand.  */
> -            int shift = clz32(abs_arg) - 9;
> -            int exp = -126 - shift + 1023;
> -            ret |= (uint64_t)exp << 52;
> -            ret |= abs_arg << (shift + 29);
> +            int lz = clz32(abs_arg);
> +            abs_arg &= ~(1 << (31 - lz)); /* [2a.] */
> +
> +            /* shift within sign32 includeing leading '1' */
> +            int shift = lz + 1 - (32 - 23);
> +            int exp = -126 + 1023 - shift; /* [2b]. */
> +            ret |= (uint64_t)exp << 52; /* [3.] */
> +            ret |= (uint64_t)abs_arg << (52 - 23 + shift); /* [4.] */
>          }
>      }
>      return ret;

-- 
David Gibson			| I'll have my music baroque, and my code
david AT gibson.dropbear.id.au	| minimalist, thank you.  NOT _the_ _other_
				| _way_ _around_!
http://www.ozlabs.org/~dgibson

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 833 bytes --]

Re: [Qemu-devel] [PATCH] powerpc: fix denorm float->double conversion

[Richard Henderson]

On 3/23/19 12:24 PM, Sergei Trofimovich wrote:
> Here denormalization conversion has a few bugs:
> - significand (abs_arg) has 32-bit unsigned wraparound in
>     ret |= abs_arg << (shift + 29);
> - significand does not drop explicit leading '1' in denorm
>   'float' when converting to normalized 'double'
> - significand had an off-by-one shift

Correct on all points.  Thanks for the test case and analysis.


> +        /*
> +         * Conversion mechanics:
> +         * float denorm (2^(-126) - biased):
> +         *    [ sign (1 bit) | exp32 (8 bits)  | sign32 (23 bits) ]
> +         *                 s                0    0001abc...def

FWIW, the overlap between "sign" and "significand" is why I prefer the term
"fraction", even though the term itself is less precise.


>          if (unlikely(abs_arg != 0)) {
>              /* Denormalized operand.  */
> -            int shift = clz32(abs_arg) - 9;
> -            int exp = -126 - shift + 1023;
> -            ret |= (uint64_t)exp << 52;
> -            ret |= abs_arg << (shift + 29);
> +            int lz = clz32(abs_arg);
> +            abs_arg &= ~(1 << (31 - lz)); /* [2a.] */
> +
> +            /* shift within sign32 includeing leading '1' */
> +            int shift = lz + 1 - (32 - 23);
> +            int exp = -126 + 1023 - shift; /* [2b]. */
> +            ret |= (uint64_t)exp << 52; /* [3.] */
> +            ret |= (uint64_t)abs_arg << (52 - 23 + shift); /* [4.] */

I think perhaps using deposit makes things clearer, since we don't have to
explicitly remove the msb in that case:

E.g.

@@ -67,10 +67,10 @@ uint64_t helper_todouble(uint32_t arg)
         ret = (uint64_t)extract32(arg, 31, 1) << 63;
         if (unlikely(abs_arg != 0)) {
             /* Denormalized operand.  */
-            int shift = clz32(abs_arg) - 9;
-            int exp = -126 - shift + 1023;
-            ret |= (uint64_t)exp << 52;
-            ret |= abs_arg << (shift + 29);
+            int msbm1 = 31 - clz32(abs_arg);
+            int exp = 1023 - 126 - (23 - msbm1);
+            ret = deposit64(ret, 52, 11, exp);
+            ret = deposit64(ret, 52 - msbm1, msbm1, abs_arg);


Thoughts?


r~

Re: [Qemu-devel] [PATCH] powerpc: fix denorm float->double conversion

[Richard Henderson]

On 3/23/19 12:24 PM, Sergei Trofimovich wrote:
> Here denormalization conversion has a few bugs:
> - significand (abs_arg) has 32-bit unsigned wraparound in
>     ret |= abs_arg << (shift + 29);
> - significand does not drop explicit leading '1' in denorm
>   'float' when converting to normalized 'double'
> - significand had an off-by-one shift

Correct on all points.  Thanks for the test case and analysis.


> +        /*
> +         * Conversion mechanics:
> +         * float denorm (2^(-126) - biased):
> +         *    [ sign (1 bit) | exp32 (8 bits)  | sign32 (23 bits) ]
> +         *                 s                0    0001abc...def

FWIW, the overlap between "sign" and "significand" is why I prefer the term
"fraction", even though the term itself is less precise.


>          if (unlikely(abs_arg != 0)) {
>              /* Denormalized operand.  */
> -            int shift = clz32(abs_arg) - 9;
> -            int exp = -126 - shift + 1023;
> -            ret |= (uint64_t)exp << 52;
> -            ret |= abs_arg << (shift + 29);
> +            int lz = clz32(abs_arg);
> +            abs_arg &= ~(1 << (31 - lz)); /* [2a.] */
> +
> +            /* shift within sign32 includeing leading '1' */
> +            int shift = lz + 1 - (32 - 23);
> +            int exp = -126 + 1023 - shift; /* [2b]. */
> +            ret |= (uint64_t)exp << 52; /* [3.] */
> +            ret |= (uint64_t)abs_arg << (52 - 23 + shift); /* [4.] */

I think perhaps using deposit makes things clearer, since we don't have to
explicitly remove the msb in that case:

E.g.

@@ -67,10 +67,10 @@ uint64_t helper_todouble(uint32_t arg)
         ret = (uint64_t)extract32(arg, 31, 1) << 63;
         if (unlikely(abs_arg != 0)) {
             /* Denormalized operand.  */
-            int shift = clz32(abs_arg) - 9;
-            int exp = -126 - shift + 1023;
-            ret |= (uint64_t)exp << 52;
-            ret |= abs_arg << (shift + 29);
+            int msbm1 = 31 - clz32(abs_arg);
+            int exp = 1023 - 126 - (23 - msbm1);
+            ret = deposit64(ret, 52, 11, exp);
+            ret = deposit64(ret, 52 - msbm1, msbm1, abs_arg);


Thoughts?


r~