[Qemu-devel] [PATCH v2 0/4] softfloat: Fix division

[Qemu-devel] [PATCH v2 0/4] softfloat: Fix division

[Richard Henderson]

Changes from v1:
  * Preserve udiv_qrnnd as a separate division primitive that
    could be used as a building block for float128 division.
  * Include asm fragments for x86_64, s390x, and ppc64.


r~


Richard Henderson (4):
  softfloat: Fix division
  softfloat: Specialize udiv_qrnnd for x86_64
  softfloat: Specialize udiv_qrnnd for s390x
  softfloat: Specialize udiv_qrnnd for ppc64

 include/fpu/softfloat-macros.h | 35 +++++++++++++++++++++++++++++++---
 fpu/softfloat.c                | 35 ++++++++++++++++++++++++++--------
 2 files changed, 59 insertions(+), 11 deletions(-)

-- 
2.17.1

[Qemu-devel] [PATCH v2 1/4] softfloat: Fix division

[Richard Henderson]

The __udiv_qrnnd primitive that we nicked from gmp requires its
inputs to be normalized.  We were not doing that.  Because the
inputs are nearly normalized already, finishing that is trivial.

Replace div128to64 with a "proper" udiv_qrnnd, so that this
remains a reusable primitive.

Fixes: cf07323d494
Fixes: https://bugs.launchpad.net/qemu/+bug/1793119
Signed-off-by: Richard Henderson <richard.henderson@linaro.org>
---
 include/fpu/softfloat-macros.h |  7 ++++---
 fpu/softfloat.c                | 35 ++++++++++++++++++++++++++--------
 2 files changed, 31 insertions(+), 11 deletions(-)

diff --git a/include/fpu/softfloat-macros.h b/include/fpu/softfloat-macros.h
index edc682139e..03312471b2 100644
--- a/include/fpu/softfloat-macros.h
+++ b/include/fpu/softfloat-macros.h
@@ -619,7 +619,8 @@ static inline uint64_t estimateDiv128To64(uint64_t a0, uint64_t a1, uint64_t b)
  *
  * Licensed under the GPLv2/LGPLv3
  */
-static inline uint64_t div128To64(uint64_t n0, uint64_t n1, uint64_t d)
+static inline uint64_t udiv_qrnnd(uint64_t *r, uint64_t n1,
+                                  uint64_t n0, uint64_t d)
 {
     uint64_t d0, d1, q0, q1, r1, r0, m;
 
@@ -658,8 +659,8 @@ static inline uint64_t div128To64(uint64_t n0, uint64_t n1, uint64_t d)
     }
     r0 -= m;
 
-    /* Return remainder in LSB */
-    return (q1 << 32) | q0 | (r0 != 0);
+    *r = r0;
+    return (q1 << 32) | q0;
 }
 
 /*----------------------------------------------------------------------------
diff --git a/fpu/softfloat.c b/fpu/softfloat.c
index a06b6ef7e4..97ef66d570 100644
--- a/fpu/softfloat.c
+++ b/fpu/softfloat.c
@@ -1112,19 +1112,38 @@ static FloatParts div_floats(FloatParts a, FloatParts b, float_status *s)
     bool sign = a.sign ^ b.sign;
 
     if (a.cls == float_class_normal && b.cls == float_class_normal) {
-        uint64_t temp_lo, temp_hi;
+        uint64_t n0, n1, q, r;
         int exp = a.exp - b.exp;
+
+        /*
+         * We want a 2*N / N-bit division to produce exactly an N-bit
+         * result, so that we do not lose any precision and so that we
+         * do not have to renormalize afterward.  If A.frac < B.frac,
+         * then division would produce an (N-1)-bit result; shift A left
+         * by one to produce the an N-bit result, and decrement the
+         * exponent to match.
+         *
+         * The udiv_qrnnd algorithm that we're using requires normalization,
+         * i.e. the msb of the denominator must be set.  Since we know that
+         * DECOMPOSED_BINARY_POINT is msb-1, the inputs must be shifted left
+         * by one (more), and the remainder must be shifted right by one.
+         */
         if (a.frac < b.frac) {
             exp -= 1;
-            shortShift128Left(0, a.frac, DECOMPOSED_BINARY_POINT + 1,
-                              &temp_hi, &temp_lo);
+            shortShift128Left(0, a.frac, DECOMPOSED_BINARY_POINT + 2, &n1, &n0);
         } else {
-            shortShift128Left(0, a.frac, DECOMPOSED_BINARY_POINT,
-                              &temp_hi, &temp_lo);
+            shortShift128Left(0, a.frac, DECOMPOSED_BINARY_POINT + 1, &n1, &n0);
         }
-        /* LSB of quot is set if inexact which roundandpack will use
-         * to set flags. Yet again we re-use a for the result */
-        a.frac = div128To64(temp_lo, temp_hi, b.frac);
+        q = udiv_qrnnd(&r, n1, n0, b.frac << 1);
+
+        /*
+         * Set lsb if there is a remainder, to set inexact.
+         * As mentioned above, to find the actual value of the remainder we
+         * would need to shift right, but (1) we are only concerned about
+         * non-zero-ness, and (2) the remainder will always be even because
+         * both inputs to the division primitive are even.
+         */
+        a.frac = q | (r != 0);
         a.sign = sign;
         a.exp = exp;
         return a;
-- 
2.17.1

[Qemu-devel] [PATCH v2 2/4] softfloat: Specialize udiv_qrnnd for x86_64

[Richard Henderson]

The ISA has a 128/64-bit division instruction.

Signed-off-by: Richard Henderson <richard.henderson@linaro.org>
---
 include/fpu/softfloat-macros.h | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/include/fpu/softfloat-macros.h b/include/fpu/softfloat-macros.h
index 03312471b2..6d58615709 100644
--- a/include/fpu/softfloat-macros.h
+++ b/include/fpu/softfloat-macros.h
@@ -622,6 +622,11 @@ static inline uint64_t estimateDiv128To64(uint64_t a0, uint64_t a1, uint64_t b)
 static inline uint64_t udiv_qrnnd(uint64_t *r, uint64_t n1,
                                   uint64_t n0, uint64_t d)
 {
+#if defined(__x86_64__)
+    uint64_t q;
+    asm("divq %4" : "=a"(q), "=d"(*r) : "0"(n0), "1"(n1), "rm"(d));
+    return q;
+#else
     uint64_t d0, d1, q0, q1, r1, r0, m;
 
     d0 = (uint32_t)d;
@@ -661,6 +666,7 @@ static inline uint64_t udiv_qrnnd(uint64_t *r, uint64_t n1,
 
     *r = r0;
     return (q1 << 32) | q0;
+#endif
 }
 
 /*----------------------------------------------------------------------------
-- 
2.17.1

[Qemu-devel] [PATCH v2 3/4] softfloat: Specialize udiv_qrnnd for s390x

[Richard Henderson]

The ISA has a 128/64-bit division instruction.

Cc: qemu-s390x@nongnu.org
Cc: Cornelia Huck <cohuck@redhat.com>
Cc: David Hildenbrand <david@redhat.com>
Signed-off-by: Richard Henderson <richard.henderson@linaro.org>
---
 include/fpu/softfloat-macros.h | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/include/fpu/softfloat-macros.h b/include/fpu/softfloat-macros.h
index 6d58615709..e702607b43 100644
--- a/include/fpu/softfloat-macros.h
+++ b/include/fpu/softfloat-macros.h
@@ -626,6 +626,12 @@ static inline uint64_t udiv_qrnnd(uint64_t *r, uint64_t n1,
     uint64_t q;
     asm("divq %4" : "=a"(q), "=d"(*r) : "0"(n0), "1"(n1), "rm"(d));
     return q;
+#elif defined(__s390x__)
+    /* Need to use a TImode type to get an even register pair for DLGR.  */
+    unsigned __int128 n = (unsigned __int128)n1 << 64 | n0;
+    asm("dlgr %0, %1" : "+r"(n) : "r"(d));
+    *r = n >> 64;
+    return n;
 #else
     uint64_t d0, d1, q0, q1, r1, r0, m;
 
-- 
2.17.1

[Qemu-devel] [PATCH v2 4/4] softfloat: Specialize udiv_qrnnd for ppc64

[Richard Henderson]

The ISA has a 128/64-bit division instruction, though it assumes the
low 64-bits of the numerator are 0, and so requires a bit more fixup
than a full 128-bit division insn.

Cc: qemu-ppc@nongnu.org
Cc: Alexander Graf <agraf@suse.de>
Cc: David Gibson <david@gibson.dropbear.id.au>
Signed-off-by: Richard Henderson <richard.henderson@linaro.org>
---
 include/fpu/softfloat-macros.h | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/include/fpu/softfloat-macros.h b/include/fpu/softfloat-macros.h
index e702607b43..001bf4f23c 100644
--- a/include/fpu/softfloat-macros.h
+++ b/include/fpu/softfloat-macros.h
@@ -632,6 +632,22 @@ static inline uint64_t udiv_qrnnd(uint64_t *r, uint64_t n1,
     asm("dlgr %0, %1" : "+r"(n) : "r"(d));
     *r = n >> 64;
     return n;
+#elif defined(_ARCH_PPC64)
+    /* From Power ISA 3.0B, programming note for divdeu.  */
+    uint64_t q1, q2, Q, r1, r2, R;
+    asm("divdeu %0,%2,%4; divdu %1,%3,%4"
+        : "=&r"(q1), "=r"(q2)
+        : "r"(n1), "r"(n0), "r"(d));
+    r1 = -(q1 * d);         /* low part of (n1<<64) - (q1 * d) */
+    r2 = n0 - (q2 * d);
+    Q = q1 + q2;
+    R = r2 + r1;
+    if (R < r2 || R >= d) { /* overflow implies R > d */
+        Q += 1;
+        R -= d;
+    }
+    *r = R;
+    return Q;
 #else
     uint64_t d0, d1, q0, q1, r1, r0, m;
 
-- 
2.17.1

Re: [Qemu-devel] [PATCH v2 4/4] softfloat: Specialize udiv_qrnnd for ppc64

[David Gibson]

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

On Wed, Oct 03, 2018 at 01:07:11PM -0500, Richard Henderson wrote:
> The ISA has a 128/64-bit division instruction, though it assumes the
> low 64-bits of the numerator are 0, and so requires a bit more fixup
> than a full 128-bit division insn.
> 
> Cc: qemu-ppc@nongnu.org
> Cc: Alexander Graf <agraf@suse.de>
> Cc: David Gibson <david@gibson.dropbear.id.au>
> Signed-off-by: Richard Henderson <richard.henderson@linaro.org>

Reviewed-by: David Gibson <david@gibson.dropbear.id.au>

> ---
>  include/fpu/softfloat-macros.h | 16 ++++++++++++++++
>  1 file changed, 16 insertions(+)
> 
> diff --git a/include/fpu/softfloat-macros.h b/include/fpu/softfloat-macros.h
> index e702607b43..001bf4f23c 100644
> --- a/include/fpu/softfloat-macros.h
> +++ b/include/fpu/softfloat-macros.h
> @@ -632,6 +632,22 @@ static inline uint64_t udiv_qrnnd(uint64_t *r, uint64_t n1,
>      asm("dlgr %0, %1" : "+r"(n) : "r"(d));
>      *r = n >> 64;
>      return n;
> +#elif defined(_ARCH_PPC64)
> +    /* From Power ISA 3.0B, programming note for divdeu.  */
> +    uint64_t q1, q2, Q, r1, r2, R;
> +    asm("divdeu %0,%2,%4; divdu %1,%3,%4"
> +        : "=&r"(q1), "=r"(q2)
> +        : "r"(n1), "r"(n0), "r"(d));
> +    r1 = -(q1 * d);         /* low part of (n1<<64) - (q1 * d) */
> +    r2 = n0 - (q2 * d);
> +    Q = q1 + q2;
> +    R = r2 + r1;
> +    if (R < r2 || R >= d) { /* overflow implies R > d */
> +        Q += 1;
> +        R -= d;
> +    }
> +    *r = R;
> +    return Q;
>  #else
>      uint64_t d0, d1, q0, q1, r1, r0, m;
>  

-- 
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 v2 0/4] softfloat: Fix division

[Alex Bennée]

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

> Changes from v1:
>   * Preserve udiv_qrnnd as a separate division primitive that
>     could be used as a building block for float128 division.
>   * Include asm fragments for x86_64, s390x, and ppc64.


It passes my fops fdiv_double test but Emilio's test is reporting:

  Errors found in f64_div, rounding near_even:
  +252.7FFFFFFFFFF80  +001.FFFFFFFFFFFFE
          => +64F.7FFFFFFFFFF82 ....x  expected +64F.7FFFFFFFFFF81 ....x
  +000.FFFFFFFFFFFFF  +3FE.FFFFFFFFFFFFF
          => +001.0000000000000 ...ux  expected +000.FFFFFFFFFFFFF ...ux
  +000.FFFFFFFFFFFFF  -3FE.FFFFFFFFFFFFF
          => -001.0000000000000 ...ux  expected -000.FFFFFFFFFFFFF ...ux
  +000.FFFFFFFFFFFFE  +3FE.FFFFFFFFFFFFF
          => +000.FFFFFFFFFFFFF ...ux  expected +000.FFFFFFFFFFFFE ...ux
  +000.FFFFFFFFFFFFE  +3FF.FFFFFFFFFFFFE
          => +000.8000000000000 ...ux  expected +000.7FFFFFFFFFFFF ...ux
  +000.FFFFFFFFFFFFE  -3FE.FFFFFFFFFFFFF
          => -000.FFFFFFFFFFFFF ...ux  expected -000.FFFFFFFFFFFFE ...ux
  +000.FFFFFFFFFFFFE  -3FF.FFFFFFFFFFFFE
          => -000.8000000000000 ...ux  expected -000.7FFFFFFFFFFFF ...ux
  -401.E7D66F623CB4C  +787.FEFF7FFFFFFFE
          => -078.E8CB4F3CAEDB8 ....x  expected -078.E8CB4F3CAEDB7 ....x
  +3FC.0000000007FEF  +198.0000002000007
          => +662.FFFFFFC00FFD1 ....x  expected +662.FFFFFFC00FFD0 ....x
  -400.FFFFF8000001F  -41C.FFFFFFBFFFFDF
          => +3E2.FFFFF84000031 ....x  expected +3E2.FFFFF84000030 ....x
  -7FE.7FFFFFFFBFFFF  +3FE.FFFFFFFFFFFFE
          => -7FE.7FFFFFFFC0001 ....x  expected -7FE.7FFFFFFFC0000 ....x
  +001.FFFFFFFFFFFFE  +3FF.FFFFFFFFFFFFF
          => +001.0000000000000 ...ux  expected +000.FFFFFFFFFFFFF ...ux
  +001.FFFFFFFFFFFFE  -3FF.FFFFFFFFFFFFF
          => -001.0000000000000 ...ux  expected -000.FFFFFFFFFFFFF ...ux
  +3CA.0000000000000  +4CF.FFFFFE000001F
          => +2F9.0000010000001 ....x  expected +2F9.0000010000000 ....x
  +3FD.0000000000000  -400.000003DFFFFFF
          => -3FB.FFFFF840001E3 ....x  expected -3FB.FFFFF840001E2 ....x
  +2C6.00000000003EF  -400.0000400400000
          => -2C4.FFFF7FF82081E ....x  expected -2C4.FFFF7FF82081D ....x
  +3FD.0000000000001  -3CA.FFFFFF0001FFF
          => -431.0000007FFF006 ....x  expected -431.0000007FFF005 ....x
  -002.FFFFFFFFFFFDE  -400.FFFFFFFFFFFFF
          => +000.FFFFFFFFFFFF0 ...ux  expected +000.FFFFFFFFFFFEF ...ux
  +3FE.0000000000001  -3FF.0000001FFFFFF
          => -3FD.FFFFFFC000005 ....x  expected -3FD.FFFFFFC000004 ....x
  +3FF.0000000000000  -47E.FFFFFC00000FF
          => -37F.000001FFFFFC1 ....x  expected -37F.000001FFFFFC0 ....x

>
> r~
>
>
> Richard Henderson (4):
>   softfloat: Fix division
>   softfloat: Specialize udiv_qrnnd for x86_64
>   softfloat: Specialize udiv_qrnnd for s390x
>   softfloat: Specialize udiv_qrnnd for ppc64
>
>  include/fpu/softfloat-macros.h | 35 +++++++++++++++++++++++++++++++---
>  fpu/softfloat.c                | 35 ++++++++++++++++++++++++++--------
>  2 files changed, 59 insertions(+), 11 deletions(-)


--
Alex Bennée

Re: [Qemu-devel] [PATCH v2 0/4] softfloat: Fix division

[Emilio G. Cota]

On Thu, Oct 04, 2018 at 10:13:55 +0100, Alex Bennée wrote:
> 
> Richard Henderson <richard.henderson@linaro.org> writes:
> 
> > Changes from v1:
> >   * Preserve udiv_qrnnd as a separate division primitive that
> >     could be used as a building block for float128 division.
> >   * Include asm fragments for x86_64, s390x, and ppc64.
> 
> It passes my fops fdiv_double test but Emilio's test is reporting:
> 
>   Errors found in f64_div, rounding near_even:
>   +252.7FFFFFFFFFF80  +001.FFFFFFFFFFFFE
>           => +64F.7FFFFFFFFFF82 ....x  expected +64F.7FFFFFFFFFF81 ....x

Did you rebuild the test program? v2 passes all f64_div tests for me.

Tested-by: Emilio G. Cota <cota@braap.org>
for patches 1 and 2.

Thanks,

		E.

Re: [Qemu-devel] [PATCH v2 0/4] softfloat: Fix division

[Alex Bennée]

Emilio G. Cota <cota@braap.org> writes:

> On Thu, Oct 04, 2018 at 10:13:55 +0100, Alex Bennée wrote:
>>
>> Richard Henderson <richard.henderson@linaro.org> writes:
>>
>> > Changes from v1:
>> >   * Preserve udiv_qrnnd as a separate division primitive that
>> >     could be used as a building block for float128 division.
>> >   * Include asm fragments for x86_64, s390x, and ppc64.
>>
>> It passes my fops fdiv_double test but Emilio's test is reporting:
>>
>>   Errors found in f64_div, rounding near_even:
>>   +252.7FFFFFFFFFF80  +001.FFFFFFFFFFFFE
>>           => +64F.7FFFFFFFFFF82 ....x  expected +64F.7FFFFFFFFFF81 ....x
>
> Did you rebuild the test program? v2 passes all f64_div tests for me.
>
> Tested-by: Emilio G. Cota <cota@braap.org>
> for patches 1 and 2.

Hmm I did:

  make clean
  make

in the tests/fp dir multiple times while going through before/after fix
scenarios. And now of course it works....

Reviewed-by: Alex Bennée <alex.bennee@linaro.org>
Tested-by: Alex Bennée <alex.bennee@linaro.org>


>
> Thanks,
>
> 		E.


--
Alex Bennée

Re: [Qemu-devel] [PATCH v2 3/4] softfloat: Specialize udiv_qrnnd for s390x

[David Hildenbrand]

On 03/10/2018 20:07, Richard Henderson wrote:
> The ISA has a 128/64-bit division instruction.
> 
> Cc: qemu-s390x@nongnu.org
> Cc: Cornelia Huck <cohuck@redhat.com>
> Cc: David Hildenbrand <david@redhat.com>
> Signed-off-by: Richard Henderson <richard.henderson@linaro.org>
> ---
>  include/fpu/softfloat-macros.h | 6 ++++++
>  1 file changed, 6 insertions(+)
> 
> diff --git a/include/fpu/softfloat-macros.h b/include/fpu/softfloat-macros.h
> index 6d58615709..e702607b43 100644
> --- a/include/fpu/softfloat-macros.h
> +++ b/include/fpu/softfloat-macros.h
> @@ -626,6 +626,12 @@ static inline uint64_t udiv_qrnnd(uint64_t *r, uint64_t n1,
>      uint64_t q;
>      asm("divq %4" : "=a"(q), "=d"(*r) : "0"(n0), "1"(n1), "rm"(d));
>      return q;
> +#elif defined(__s390x__)
> +    /* Need to use a TImode type to get an even register pair for DLGR.  */
> +    unsigned __int128 n = (unsigned __int128)n1 << 64 | n0;
> +    asm("dlgr %0, %1" : "+r"(n) : "r"(d));
> +    *r = n >> 64;
> +    return n;
>  #else
>      uint64_t d0, d1, q0, q1, r1, r0, m;
>  
> 

Reviewed-by: David Hildenbrand <david@redhat.com>

-- 

Thanks,

David / dhildenb