This is a little bit of a departure from softfloat's original approach
as we skip the estimate step in favour of a straight iteration.

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

---
v3
  - added to series
  - fixed renames of structs
v4
  - fix up comments
  - use is_nan
  - use return_nan instead of pick_nan(a,a)
---
 fpu/softfloat.c         | 201 ++++++++++++++++++++++--------------------------
 include/fpu/softfloat.h |   1 +
 2 files changed, 91 insertions(+), 111 deletions(-)

diff --git a/fpu/softfloat.c b/fpu/softfloat.c
index 8fc1c2a8d9..80301d8e04 100644
--- a/fpu/softfloat.c
+++ b/fpu/softfloat.c
@@ -1897,6 +1897,96 @@ float64 float64_scalbn(float64 a, int n, float_status *status)
     return float64_round_pack_canonical(pr, status);
 }
 
+/*
+ * Square Root
+ *
+ * The old softfloat code did an approximation step before zeroing in
+ * on the final result. However for simpleness we just compute the
+ * square root by iterating down from the implicit bit to enough extra
+ * bits to ensure we get a correctly rounded result.
+ */
+
+static FloatParts sqrt_float(FloatParts a, float_status *s, const FloatFmt *p)
+{
+    uint64_t a_frac, r_frac, s_frac;
+    int bit, last_bit;
+
+    if (is_nan(a.cls)) {
+        return return_nan(a, s);
+    }
+    if (a.cls == float_class_zero) {
+        return a;  /* sqrt(+-0) = +-0 */
+    }
+    if (a.sign) {
+        s->float_exception_flags |= float_flag_invalid;
+        a.cls = float_class_dnan;
+        return a;
+    }
+    if (a.cls == float_class_inf) {
+        return a;  /* sqrt(+inf) = +inf */
+    }
+
+    assert(a.cls == float_class_normal);
+
+    /* We need two overflow bits at the top.  Adding room for that is
+       a right shift.  If the exponent is odd, we can discard the low
+       bit by multiplying the fraction by 2; that's a left shift.
+       Combine those and we shift right if the exponent is even.  */
+    a_frac = a.frac;
+    if (!(a.exp & 1)) {
+        a_frac >>= 1;
+    }
+    a.exp >>= 1;
+
+    /* Bit-by-bit computation of sqrt.  */
+    r_frac = 0;
+    s_frac = 0;
+
+    /* Iterate from implicit bit down to the 3 extra bits to compute a
+     * properly rounded result.  Remember we've inserted one more bit
+     * at the top, so these positions are one less.  */
+    bit = DECOMPOSED_BINARY_POINT - 1;
+    last_bit = MAX(p->frac_shift - 4, 0);
+    do {
+        uint64_t q = 1ULL << bit;
+        uint64_t t_frac = s_frac + q;
+        if (t_frac <= a_frac) {
+            s_frac = t_frac + q;
+            a_frac -= t_frac;
+            r_frac += q;
+        }
+        a_frac <<= 1;
+    } while (--bit >= last_bit);
+
+    /* Undo the right shift done above.  If there is any remaining
+       fraction, the result is inexact.  Set the sticky bit.  */
+    a.frac = (r_frac << 1) + (a_frac != 0);
+
+    return a;
+}
+
+float16 float16_sqrt(float16 a, float_status *status)
+{
+    FloatParts pa = float16_unpack_canonical(a, status);
+    FloatParts pr = sqrt_float(pa, status, &float16_params);
+    return float16_round_pack_canonical(pr, status);
+}
+
+float32 float32_sqrt(float32 a, float_status *status)
+{
+    FloatParts pa = float32_unpack_canonical(a, status);
+    FloatParts pr = sqrt_float(pa, status, &float32_params);
+    return float32_round_pack_canonical(pr, status);
+}
+
+float64 float64_sqrt(float64 a, float_status *status)
+{
+    FloatParts pa = float64_unpack_canonical(a, status);
+    FloatParts pr = sqrt_float(pa, status, &float64_params);
+    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
@@ -3304,62 +3394,6 @@ float32 float32_rem(float32 a, float32 b, float_status *status)
 }
 
 
-/*----------------------------------------------------------------------------
-| Returns the square root of the single-precision floating-point value `a'.
-| The operation is performed according to the IEC/IEEE Standard for Binary
-| Floating-Point Arithmetic.
-*----------------------------------------------------------------------------*/
-
-float32 float32_sqrt(float32 a, float_status *status)
-{
-    flag aSign;
-    int aExp, zExp;
-    uint32_t aSig, zSig;
-    uint64_t rem, term;
-    a = float32_squash_input_denormal(a, status);
-
-    aSig = extractFloat32Frac( a );
-    aExp = extractFloat32Exp( a );
-    aSign = extractFloat32Sign( a );
-    if ( aExp == 0xFF ) {
-        if (aSig) {
-            return propagateFloat32NaN(a, float32_zero, status);
-        }
-        if ( ! aSign ) return a;
-        float_raise(float_flag_invalid, status);
-        return float32_default_nan(status);
-    }
-    if ( aSign ) {
-        if ( ( aExp | aSig ) == 0 ) return a;
-        float_raise(float_flag_invalid, status);
-        return float32_default_nan(status);
-    }
-    if ( aExp == 0 ) {
-        if ( aSig == 0 ) return float32_zero;
-        normalizeFloat32Subnormal( aSig, &aExp, &aSig );
-    }
-    zExp = ( ( aExp - 0x7F )>>1 ) + 0x7E;
-    aSig = ( aSig | 0x00800000 )<<8;
-    zSig = estimateSqrt32( aExp, aSig ) + 2;
-    if ( ( zSig & 0x7F ) <= 5 ) {
-        if ( zSig < 2 ) {
-            zSig = 0x7FFFFFFF;
-            goto roundAndPack;
-        }
-        aSig >>= aExp & 1;
-        term = ( (uint64_t) zSig ) * zSig;
-        rem = ( ( (uint64_t) aSig )<<32 ) - term;
-        while ( (int64_t) rem < 0 ) {
-            --zSig;
-            rem += ( ( (uint64_t) zSig )<<1 ) | 1;
-        }
-        zSig |= ( rem != 0 );
-    }
-    shift32RightJamming( zSig, 1, &zSig );
- roundAndPack:
-    return roundAndPackFloat32(0, zExp, zSig, status);
-
-}
 
 /*----------------------------------------------------------------------------
 | Returns the binary exponential of the single-precision floating-point value
@@ -4203,61 +4237,6 @@ float64 float64_rem(float64 a, float64 b, float_status *status)
 
 }
 
-
-/*----------------------------------------------------------------------------
-| Returns the square root of the double-precision floating-point value `a'.
-| The operation is performed according to the IEC/IEEE Standard for Binary
-| Floating-Point Arithmetic.
-*----------------------------------------------------------------------------*/
-
-float64 float64_sqrt(float64 a, float_status *status)
-{
-    flag aSign;
-    int aExp, zExp;
-    uint64_t aSig, zSig, doubleZSig;
-    uint64_t rem0, rem1, term0, term1;
-    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);
-        }
-        if ( ! aSign ) return a;
-        float_raise(float_flag_invalid, status);
-        return float64_default_nan(status);
-    }
-    if ( aSign ) {
-        if ( ( aExp | aSig ) == 0 ) return a;
-        float_raise(float_flag_invalid, status);
-        return float64_default_nan(status);
-    }
-    if ( aExp == 0 ) {
-        if ( aSig == 0 ) return float64_zero;
-        normalizeFloat64Subnormal( aSig, &aExp, &aSig );
-    }
-    zExp = ( ( aExp - 0x3FF )>>1 ) + 0x3FE;
-    aSig |= LIT64( 0x0010000000000000 );
-    zSig = estimateSqrt32( aExp, aSig>>21 );
-    aSig <<= 9 - ( aExp & 1 );
-    zSig = estimateDiv128To64( aSig, 0, zSig<<32 ) + ( zSig<<30 );
-    if ( ( zSig & 0x1FF ) <= 5 ) {
-        doubleZSig = zSig<<1;
-        mul64To128( zSig, zSig, &term0, &term1 );
-        sub128( aSig, 0, term0, term1, &rem0, &rem1 );
-        while ( (int64_t) rem0 < 0 ) {
-            --zSig;
-            doubleZSig -= 2;
-            add128( rem0, rem1, zSig>>63, doubleZSig | 1, &rem0, &rem1 );
-        }
-        zSig |= ( ( rem0 | rem1 ) != 0 );
-    }
-    return roundAndPackFloat64(0, zExp, zSig, status);
-
-}
-
 /*----------------------------------------------------------------------------
 | Returns the binary log of the double-precision floating-point value `a'.
 | The operation is performed according to the IEC/IEEE Standard for Binary
diff --git a/include/fpu/softfloat.h b/include/fpu/softfloat.h
index cebe37b716..9b7b5e34e2 100644
--- a/include/fpu/softfloat.h
+++ b/include/fpu/softfloat.h
@@ -251,6 +251,7 @@ float16 float16_minnum(float16, float16, float_status *status);
 float16 float16_maxnum(float16, float16, float_status *status);
 float16 float16_minnummag(float16, float16, float_status *status);
 float16 float16_maxnummag(float16, float16, float_status *status);
+float16 float16_sqrt(float16, float_status *status);
 int float16_compare(float16, float16, float_status *status);
 int float16_compare_quiet(float16, float16, float_status *status);
 
-- 
2.15.1

On 6 February 2018 at 16:48, Alex Bennée <alex.bennee@linaro.org> wrote:
> This is a little bit of a departure from softfloat's original approach
> as we skip the estimate step in favour of a straight iteration.
>
> Suggested-by: Richard Henderson <richard.henderson@linaro.org>
> Signed-off-by: Alex Bennée <alex.bennee@linaro.org>
>
> ---
> v3
>   - added to series
>   - fixed renames of structs
> v4
>   - fix up comments
>   - use is_nan
>   - use return_nan instead of pick_nan(a,a)
> ---
>  fpu/softfloat.c         | 201 ++++++++++++++++++++++--------------------------
>  include/fpu/softfloat.h |   1 +
>  2 files changed, 91 insertions(+), 111 deletions(-)
>
> diff --git a/fpu/softfloat.c b/fpu/softfloat.c
> index 8fc1c2a8d9..80301d8e04 100644
> --- a/fpu/softfloat.c
> +++ b/fpu/softfloat.c
> @@ -1897,6 +1897,96 @@ float64 float64_scalbn(float64 a, int n, float_status *status)
>      return float64_round_pack_canonical(pr, status);
>  }
>
> +/*
> + * Square Root
> + *
> + * The old softfloat code did an approximation step before zeroing in
> + * on the final result. However for simpleness we just compute the
> + * square root by iterating down from the implicit bit to enough extra
> + * bits to ensure we get a correctly rounded result.

So is this new approach slower?

> + */
> +
> +static FloatParts sqrt_float(FloatParts a, float_status *s, const FloatFmt *p)
> +{
> +    uint64_t a_frac, r_frac, s_frac;
> +    int bit, last_bit;
> +
> +    if (is_nan(a.cls)) {
> +        return return_nan(a, s);
> +    }
> +    if (a.cls == float_class_zero) {
> +        return a;  /* sqrt(+-0) = +-0 */
> +    }
> +    if (a.sign) {
> +        s->float_exception_flags |= float_flag_invalid;
> +        a.cls = float_class_dnan;
> +        return a;
> +    }
> +    if (a.cls == float_class_inf) {
> +        return a;  /* sqrt(+inf) = +inf */
> +    }
> +
> +    assert(a.cls == float_class_normal);
> +
> +    /* We need two overflow bits at the top.  Adding room for that is
> +       a right shift.  If the exponent is odd, we can discard the low
> +       bit by multiplying the fraction by 2; that's a left shift.
> +       Combine those and we shift right if the exponent is even.  */
> +    a_frac = a.frac;
> +    if (!(a.exp & 1)) {
> +        a_frac >>= 1;
> +    }
> +    a.exp >>= 1;

Comment says "shift right if the exponent is even", but code
says "shift right by 1 if exponent is odd, by 2 if exponent is even".

> +
> +    /* Bit-by-bit computation of sqrt.  */
> +    r_frac = 0;
> +    s_frac = 0;
> +
> +    /* Iterate from implicit bit down to the 3 extra bits to compute a
> +     * properly rounded result.  Remember we've inserted one more bit
> +     * at the top, so these positions are one less.  */

Some consistency in multiline comment formatting would be nice.
The top-of-function header and these two in the body of
the function have 3 different styles between them...

Otherwise
Reviewed-by: Peter Maydell <peter.maydell@linaro.org>

thanks
-- PMM

On 02/13/2018 07:50 AM, Peter Maydell wrote:
>> +    /* We need two overflow bits at the top.  Adding room for that is
>> +       a right shift.  If the exponent is odd, we can discard the low
>> +       bit by multiplying the fraction by 2; that's a left shift.
>> +       Combine those and we shift right if the exponent is even.  */
>> +    a_frac = a.frac;
>> +    if (!(a.exp & 1)) {
>> +        a_frac >>= 1;
>> +    }
>> +    a.exp >>= 1;
> Comment says "shift right if the exponent is even", but code
> says "shift right by 1 if exponent is odd, by 2 if exponent is even".
> 

The last line is dividing the exponent by 2, not shifting the fraction.


r~

On 13 February 2018 at 16:23, Richard Henderson
<richard.henderson@linaro.org> wrote:
> On 02/13/2018 07:50 AM, Peter Maydell wrote:
>>> +    /* We need two overflow bits at the top.  Adding room for that is
>>> +       a right shift.  If the exponent is odd, we can discard the low
>>> +       bit by multiplying the fraction by 2; that's a left shift.
>>> +       Combine those and we shift right if the exponent is even.  */
>>> +    a_frac = a.frac;
>>> +    if (!(a.exp & 1)) {
>>> +        a_frac >>= 1;
>>> +    }
>>> +    a.exp >>= 1;
>> Comment says "shift right if the exponent is even", but code
>> says "shift right by 1 if exponent is odd, by 2 if exponent is even".
>>
>
> The last line is dividing the exponent by 2, not shifting the fraction.

Doh, so it is.

-- PMM

This is an attempt to save some of the cost of sqrt by using the
inbuilt support of the host hardware. The idea is assuming we start
with a valid input we can use the hardware. If any tininess issues
occur this will trip and FPU exception where:

  - we turn off cpu->use_host_fpu
  - mask the FPU exceptions
  - return to what we were doing

Once we return we should pick up the fact that there was something
weird about the operation and fall-back to the pure software
implementation.

You could imagine this being extended for code generation but instead
of returning to the code we could exit and re-generate the TB but this
time with pure software helpers rather than any support from the
hardware.

This is a sort of fix-it-up after the fact approach because reading
the FP state is an expensive operation for everything so let's only
worry about exceptions when they trip...

Signed-off-by: Alex Bennée <alex.bennee@linaro.org>
---
 cpus.c                        | 28 ++++++++++++++++++++++++++++
 fpu/softfloat.c               | 40 +++++++++++++++++++++++++++++++++++-----
 include/fpu/softfloat-types.h |  2 ++
 include/fpu/softfloat.h       |  4 ++++
 include/qom/cpu.h             |  1 +
 linux-user/main.c             |  8 ++++++++
 linux-user/signal.c           | 16 ++++++++++++++++
 target/arm/cpu.c              |  4 ++++
 8 files changed, 98 insertions(+), 5 deletions(-)

diff --git a/cpus.c b/cpus.c
index f298b659f4..e435f6737b 100644
--- a/cpus.c
+++ b/cpus.c
@@ -23,6 +23,7 @@
  */
 
 #include "qemu/osdep.h"
+#include <fenv.h>
 #include "qemu/config-file.h"
 #include "cpu.h"
 #include "monitor/monitor.h"
@@ -1078,10 +1079,36 @@ static void qemu_init_sigbus(void)
 
     prctl(PR_MCE_KILL, PR_MCE_KILL_SET, PR_MCE_KILL_EARLY, 0, 0);
 }
+
+static void sigfpu_handler(int n, siginfo_t *siginfo, void *ctx)
+{
+    fprintf(stderr, "%s: got %d, %p/%p\n", __func__, n, siginfo, ctx);
+
+    /* Called asynchronously in VCPU thread.  */
+    g_assert(current_cpu);
+}
+
+static void qemu_init_sigfpu(void)
+{
+    struct sigaction action;
+
+    memset(&action, 0, sizeof(action));
+    action.sa_flags = SA_SIGINFO;
+    action.sa_sigaction = sigfpu_handler;
+    sigaction(SIGBUS, &action, NULL);
+
+    feenableexcept(FE_INVALID   |
+                   FE_OVERFLOW  |
+                   FE_UNDERFLOW |
+                   FE_INEXACT);
+}
 #else /* !CONFIG_LINUX */
 static void qemu_init_sigbus(void)
 {
 }
+static void qemu_init_sigfpu(void)
+{
+}
 #endif /* !CONFIG_LINUX */
 
 static QemuMutex qemu_global_mutex;
@@ -1827,6 +1854,7 @@ static void qemu_tcg_init_vcpu(CPUState *cpu)
     if (!tcg_region_inited) {
         tcg_region_inited = 1;
         tcg_region_init();
+        qemu_init_sigfpu();
     }
 
     if (qemu_tcg_mttcg_enabled() || !single_tcg_cpu_thread) {
diff --git a/fpu/softfloat.c b/fpu/softfloat.c
index e7fb0d357a..ec9355af7a 100644
--- a/fpu/softfloat.c
+++ b/fpu/softfloat.c
@@ -1905,10 +1905,12 @@ float64 float64_scalbn(float64 a, int n, float_status *status)
  * bits to ensure we get a correctly rounded result.
  *
  * This does mean however the calculation is slower than before,
- * especially for 64 bit floats.
+ * especially for 64 bit floats. However the caller can only do checks
+ * if they actually want to off-load to the library.
  */
 
-static FloatParts sqrt_float(FloatParts a, float_status *s, const FloatFmt *p)
+static FloatParts sqrt_float(FloatParts a, float_status *s,
+                             const FloatFmt *p, bool check_only)
 {
     uint64_t a_frac, r_frac, s_frac;
     int bit, last_bit;
@@ -1928,6 +1930,10 @@ static FloatParts sqrt_float(FloatParts a, float_status *s, const FloatFmt *p)
         return a;  /* sqrt(+inf) = +inf */
     }
 
+    if (check_only) {
+        return a;
+    }
+
     assert(a.cls == float_class_normal);
 
     /* We need two overflow bits at the top. Adding room for that is a
@@ -1973,21 +1979,45 @@ static FloatParts sqrt_float(FloatParts a, float_status *s, const FloatFmt *p)
 float16 __attribute__((flatten)) float16_sqrt(float16 a, float_status *status)
 {
     FloatParts pa = float16_unpack_canonical(a, status);
-    FloatParts pr = sqrt_float(pa, status, &float16_params);
+    FloatParts pr = sqrt_float(pa, status, &float16_params, false);
     return float16_round_pack_canonical(pr, status);
 }
 
 float32 __attribute__((flatten)) float32_sqrt(float32 a, float_status *status)
 {
     FloatParts pa = float32_unpack_canonical(a, status);
-    FloatParts pr = sqrt_float(pa, status, &float32_params);
+    FloatParts pr;
+
+    if (status->use_host_fpu && *status->use_host_fpu) {
+        pr = sqrt_float(pa, status, &float32_params, true);
+        if (pr.cls == float_class_normal) {
+            float32 r = __builtin_sqrt(a);
+            if (*status->use_host_fpu) {
+                return r;
+            }
+        }
+    }
+
+    pr = sqrt_float(pa, status, &float32_params, false);
     return float32_round_pack_canonical(pr, status);
 }
 
 float64 __attribute__((flatten)) float64_sqrt(float64 a, float_status *status)
 {
     FloatParts pa = float64_unpack_canonical(a, status);
-    FloatParts pr = sqrt_float(pa, status, &float64_params);
+    FloatParts pr = sqrt_float(pa, status, &float64_params, true);
+
+    if (status->use_host_fpu && *status->use_host_fpu) {
+        pr = sqrt_float(pa, status, &float64_params, true);
+        if (pr.cls == float_class_normal) {
+            float64 r = __builtin_sqrt(a);
+            if (*status->use_host_fpu) {
+                return r;
+            }
+        }
+    }
+
+    pr = sqrt_float(pa, status, &float64_params, false);
     return float64_round_pack_canonical(pr, status);
 }
 
diff --git a/include/fpu/softfloat-types.h b/include/fpu/softfloat-types.h
index 4e378cb612..4c32e56cad 100644
--- a/include/fpu/softfloat-types.h
+++ b/include/fpu/softfloat-types.h
@@ -174,6 +174,8 @@ typedef struct float_status {
     flag flush_inputs_to_zero;
     flag default_nan_mode;
     flag snan_bit_is_one;
+    /* can we use the host_fpu for some things? */
+    bool *use_host_fpu;
 } float_status;
 
 #endif /* SOFTFLOAT_TYPES_H */
diff --git a/include/fpu/softfloat.h b/include/fpu/softfloat.h
index 9b7b5e34e2..f7ee0232a2 100644
--- a/include/fpu/softfloat.h
+++ b/include/fpu/softfloat.h
@@ -157,6 +157,10 @@ static inline flag get_default_nan_mode(float_status *status)
 {
     return status->default_nan_mode;
 }
+static inline void enable_host_fpu(bool *host_fpu_flag, float_status *status)
+{
+    status->use_host_fpu = host_fpu_flag;
+}
 
 /*----------------------------------------------------------------------------
 | Routine to raise any or all of the software IEC/IEEE floating-point
diff --git a/include/qom/cpu.h b/include/qom/cpu.h
index aff88fa16f..337ebef8b6 100644
--- a/include/qom/cpu.h
+++ b/include/qom/cpu.h
@@ -396,6 +396,7 @@ struct CPUState {
     uint32_t halted;
     uint32_t can_do_io;
     int32_t exception_index;
+    bool use_host_fpu;
 
     /* shared by kvm, hax and hvf */
     bool vcpu_dirty;
diff --git a/linux-user/main.c b/linux-user/main.c
index 7de0e02487..36b6be3b2b 100644
--- a/linux-user/main.c
+++ b/linux-user/main.c
@@ -20,6 +20,7 @@
 #include "qemu-version.h"
 #include <sys/syscall.h>
 #include <sys/resource.h>
+#include <fenv.h>
 
 #include "qapi/error.h"
 #include "qemu.h"
@@ -4927,6 +4928,13 @@ int main(int argc, char **argv, char **envp)
         }
         gdb_handlesig(cpu, 0);
     }
+
+    feenableexcept(FE_INVALID   |
+                   FE_OVERFLOW  |
+                   FE_UNDERFLOW |
+                   FE_INEXACT);
+    cpu->use_host_fpu = true;
+
     cpu_loop(env);
     /* never exits */
     return 0;
diff --git a/linux-user/signal.c b/linux-user/signal.c
index 9a380b9e31..0773d3ef18 100644
--- a/linux-user/signal.c
+++ b/linux-user/signal.c
@@ -20,6 +20,7 @@
 #include "qemu/bitops.h"
 #include <sys/ucontext.h>
 #include <sys/resource.h>
+#include <fenv.h>
 
 #include "qemu.h"
 #include "qemu-common.h"
@@ -639,6 +640,21 @@ static void host_signal_handler(int host_signum, siginfo_t *info,
     ucontext_t *uc = puc;
     struct emulated_sigtable *k;
 
+    /* Catch any FPU exceptions we might get from having tried to use
+     * the host FPU to speed up some calculations
+     */
+    if (host_signum == SIGFPE && cpu->use_host_fpu) {
+        cpu->use_host_fpu = false;
+        /* sadly this gets lost on the context switch when we return */
+        fedisableexcept(FE_INVALID   |
+                        FE_OVERFLOW  |
+                        FE_UNDERFLOW |
+                        FE_INEXACT);
+        /* sigaddset(&uc->uc_sigmask, SIGFPE); */
+        uc->__fpregs_mem.mxcsr |= 0x1f80;
+        return;
+    }
+
     /* the CPU emulator uses some host signals to detect exceptions,
        we forward to it some signals */
     if ((host_signum == SIGSEGV || host_signum == SIGBUS)
diff --git a/target/arm/cpu.c b/target/arm/cpu.c
index 1b3ae62db6..67dce53a68 100644
--- a/target/arm/cpu.c
+++ b/target/arm/cpu.c
@@ -306,6 +306,10 @@ static void arm_cpu_reset(CPUState *s)
                               &env->vfp.fp_status);
     set_float_detect_tininess(float_tininess_before_rounding,
                               &env->vfp.standard_fp_status);
+
+    enable_host_fpu(&s->use_host_fpu, &env->vfp.fp_status);
+    enable_host_fpu(&s->use_host_fpu, &env->vfp.standard_fp_status);
+
 #ifndef CONFIG_USER_ONLY
     if (kvm_enabled()) {
         kvm_arm_reset_vcpu(cpu);
-- 
2.15.1

Alex Bennée <alex.bennee@linaro.org> writes:

> This is an attempt to save some of the cost of sqrt by using the
> inbuilt support of the host hardware. The idea is assuming we start
> with a valid input we can use the hardware. If any tininess issues
> occur this will trip and FPU exception where:
>
>   - we turn off cpu->use_host_fpu
>   - mask the FPU exceptions
>   - return to what we were doing
>
> Once we return we should pick up the fact that there was something
> weird about the operation and fall-back to the pure software
> implementation.
>
> You could imagine this being extended for code generation but instead
> of returning to the code we could exit and re-generate the TB but this
> time with pure software helpers rather than any support from the
> hardware.
>
> This is a sort of fix-it-up after the fact approach because reading
> the FP state is an expensive operation for everything so let's only
> worry about exceptions when they trip...
>
<snip>
> --- a/linux-user/signal.c
> +++ b/linux-user/signal.c
> @@ -20,6 +20,7 @@
>  #include "qemu/bitops.h"
>  #include <sys/ucontext.h>
>  #include <sys/resource.h>
> +#include <fenv.h>
>
>  #include "qemu.h"
>  #include "qemu-common.h"
> @@ -639,6 +640,21 @@ static void host_signal_handler(int host_signum, siginfo_t *info,
>      ucontext_t *uc = puc;
>      struct emulated_sigtable *k;
>
> +    /* Catch any FPU exceptions we might get from having tried to use
> +     * the host FPU to speed up some calculations
> +     */
> +    if (host_signum == SIGFPE && cpu->use_host_fpu) {
> +        cpu->use_host_fpu = false;
> +        /* sadly this gets lost on the context switch when we return */
> +        fedisableexcept(FE_INVALID   |
> +                        FE_OVERFLOW  |
> +                        FE_UNDERFLOW |
> +                        FE_INEXACT);
> +        /* sigaddset(&uc->uc_sigmask, SIGFPE); */
> +        uc->__fpregs_mem.mxcsr |= 0x1f80;

This is a bug, the correct place to reset mxcsr for the return is:

        (uc->uc_mcontext.fpregs)->mxcsr |= 0x1f80;

--
Alex Bennée

On Tue, Feb 20, 2018 at 21:01:37 +0000, Alex Bennée wrote:
> This is an attempt to save some of the cost of sqrt by using the
> inbuilt support of the host hardware. The idea is assuming we start
> with a valid input we can use the hardware. If any tininess issues
> occur this will trip and FPU exception where:
> 
>   - we turn off cpu->use_host_fpu
>   - mask the FPU exceptions
>   - return to what we were doing

Please see this thread for an alternative approach:

  [Qemu-devel] [PATCH v1 00/14] fp-test + hostfloat
  https://lists.gnu.org/archive/html/qemu-devel/2018-03/msg05908.html

		Emilio

On 02/06/2018 08:48 AM, Alex Bennée wrote:
> This is a little bit of a departure from softfloat's original approach
> as we skip the estimate step in favour of a straight iteration.
> 
> Suggested-by: Richard Henderson <richard.henderson@linaro.org>
> Signed-off-by: Alex Bennée <alex.bennee@linaro.org>
> 
> ---
> v3
>   - added to series
>   - fixed renames of structs
> v4
>   - fix up comments
>   - use is_nan
>   - use return_nan instead of pick_nan(a,a)
> ---
>  fpu/softfloat.c         | 201 ++++++++++++++++++++++--------------------------
>  include/fpu/softfloat.h |   1 +
>  2 files changed, 91 insertions(+), 111 deletions(-)

Did I not give you an r-b for this before?  Anyway,

Reviewed-by: Richard Henderson <richard.henderson@linaro.org>


r~