target-arm: A64: Implement AdvSIMD reciprocal estimate insns URECPE, FRECPE

Implement URECPE and FRECPE instructions in both scalar and vector forms.
The actual reciprocal estimate function is shared with the A32/T32 Neon
code. However in A64 we aren't using the Neon "standard FPSCR value"
so extra checks are necessary to handle non-squashed denormal inputs
which can never happen for A32/T32. Calling conventions for the helpers
are thus modified to pass the fpst directly; we mark the helpers as
TCG_CALL_NO_RWG since we're changing the declarations anyway.

Signed-off-by: Alex Bennée <alex.bennee@linaro.org>
Signed-off-by: Peter Maydell <peter.maydell@linaro.org>
Reviewed-by: Richard Henderson <rth@twiddle.net>
Message-id: 1394822294-14837-21-git-send-email-peter.maydell@linaro.org

1 files changed, 163 insertions, 35 deletions

diff --git a/target-arm/helper.c b/target-arm/helper.c
index 5080372bad..535fc8f35e 100644
--- a/target-arm/helper.c
+++ b/target-arm/helper.c
@@ -4520,16 +4520,21 @@ float32 HELPER(rsqrts_f32)(float32 a, float32 b, CPUARMState *env)
  * int->float conversions at run-time.  */
 #define float64_256 make_float64(0x4070000000000000LL)
 #define float64_512 make_float64(0x4080000000000000LL)
+#define float32_maxnorm make_float32(0x7f7fffff)
+#define float64_maxnorm make_float64(0x7fefffffffffffffLL)
 
-/* The algorithm that must be used to calculate the estimate
- * is specified by the ARM ARM.
+/* Reciprocal functions
+ *
+ * The algorithm that must be used to calculate the estimate
+ * is specified by the ARM ARM, see FPRecipEstimate()
  */
-static float64 recip_estimate(float64 a, CPUARMState *env)
+
+static float64 recip_estimate(float64 a, float_status *real_fp_status)
 {
     /* These calculations mustn't set any fp exception flags,
      * so we use a local copy of the fp_status.
      */
-    float_status dummy_status = env->vfp.standard_fp_status;
+    float_status dummy_status = *real_fp_status;
     float_status *s = &dummy_status;
     /* q = (int)(a * 512.0) */
     float64 q = float64_mul(float64_512, a, s);
@@ -4550,45 +4555,167 @@ static float64 recip_estimate(float64 a, CPUARMState *env)
     return float64_div(int64_to_float64(q_int, s), float64_256, s);
 }
 
-float32 HELPER(recpe_f32)(float32 a, CPUARMState *env)
+/* Common wrapper to call recip_estimate */
+static float64 call_recip_estimate(float64 num, int off, float_status *fpst)
 {
-    float_status *s = &env->vfp.standard_fp_status;
-    float64 f64;
-    uint32_t val32 = float32_val(a);
+    uint64_t val64 = float64_val(num);
+    uint64_t frac = extract64(val64, 0, 52);
+    int64_t exp = extract64(val64, 52, 11);
+    uint64_t sbit;
+    float64 scaled, estimate;
 
-    int result_exp;
-    int a_exp = (val32  & 0x7f800000) >> 23;
-    int sign = val32 & 0x80000000;
+    /* Generate the scaled number for the estimate function */
+    if (exp == 0) {
+        if (extract64(frac, 51, 1) == 0) {
+            exp = -1;
+            frac = extract64(frac, 0, 50) << 2;
+        } else {
+            frac = extract64(frac, 0, 51) << 1;
+        }
+    }
 
-    if (float32_is_any_nan(a)) {
-        if (float32_is_signaling_nan(a)) {
-            float_raise(float_flag_invalid, s);
+    /* scaled = '0' : '01111111110' : fraction<51:44> : Zeros(44); */
+    scaled = make_float64((0x3feULL << 52)
+                          | extract64(frac, 44, 8) << 44);
+
+    estimate = recip_estimate(scaled, fpst);
+
+    /* Build new result */
+    val64 = float64_val(estimate);
+    sbit = 0x8000000000000000ULL & val64;
+    exp = off - exp;
+    frac = extract64(val64, 0, 52);
+
+    if (exp == 0) {
+        frac = 1ULL << 51 | extract64(frac, 1, 51);
+    } else if (exp == -1) {
+        frac = 1ULL << 50 | extract64(frac, 2, 50);
+        exp = 0;
+    }
+
+    return make_float64(sbit | (exp << 52) | frac);
+}
+
+static bool round_to_inf(float_status *fpst, bool sign_bit)
+{
+    switch (fpst->float_rounding_mode) {
+    case float_round_nearest_even: /* Round to Nearest */
+        return true;
+    case float_round_up: /* Round to +Inf */
+        return !sign_bit;
+    case float_round_down: /* Round to -Inf */
+        return sign_bit;
+    case float_round_to_zero: /* Round to Zero */
+        return false;
+    }
+
+    g_assert_not_reached();
+}
+
+float32 HELPER(recpe_f32)(float32 input, void *fpstp)
+{
+    float_status *fpst = fpstp;
+    float32 f32 = float32_squash_input_denormal(input, fpst);
+    uint32_t f32_val = float32_val(f32);
+    uint32_t f32_sbit = 0x80000000ULL & f32_val;
+    int32_t f32_exp = extract32(f32_val, 23, 8);
+    uint32_t f32_frac = extract32(f32_val, 0, 23);
+    float64 f64, r64;
+    uint64_t r64_val;
+    int64_t r64_exp;
+    uint64_t r64_frac;
+
+    if (float32_is_any_nan(f32)) {
+        float32 nan = f32;
+        if (float32_is_signaling_nan(f32)) {
+            float_raise(float_flag_invalid, fpst);
+            nan = float32_maybe_silence_nan(f32);
         }
-        return float32_default_nan;
-    } else if (float32_is_infinity(a)) {
-        return float32_set_sign(float32_zero, float32_is_neg(a));
-    } else if (float32_is_zero_or_denormal(a)) {
-        if (!float32_is_zero(a)) {
-            float_raise(float_flag_input_denormal, s);
+        if (fpst->default_nan_mode) {
+            nan =  float32_default_nan;
         }
-        float_raise(float_flag_divbyzero, s);
-        return float32_set_sign(float32_infinity, float32_is_neg(a));
-    } else if (a_exp >= 253) {
-        float_raise(float_flag_underflow, s);
-        return float32_set_sign(float32_zero, float32_is_neg(a));
+        return nan;
+    } else if (float32_is_infinity(f32)) {
+        return float32_set_sign(float32_zero, float32_is_neg(f32));
+    } else if (float32_is_zero(f32)) {
+        float_raise(float_flag_divbyzero, fpst);
+        return float32_set_sign(float32_infinity, float32_is_neg(f32));
+    } else if ((f32_val & ~(1ULL << 31)) < (1ULL << 21)) {
+        /* Abs(value) < 2.0^-128 */
+        float_raise(float_flag_overflow | float_flag_inexact, fpst);
+        if (round_to_inf(fpst, f32_sbit)) {
+            return float32_set_sign(float32_infinity, float32_is_neg(f32));
+        } else {
+            return float32_set_sign(float32_maxnorm, float32_is_neg(f32));
+        }
+    } else if (f32_exp >= 253 && fpst->flush_to_zero) {
+        float_raise(float_flag_underflow, fpst);
+        return float32_set_sign(float32_zero, float32_is_neg(f32));
     }
 
-    f64 = make_float64((0x3feULL << 52)
-                       | ((int64_t)(val32 & 0x7fffff) << 29));
 
-    result_exp = 253 - a_exp;
+    f64 = make_float64(((int64_t)(f32_exp) << 52) | (int64_t)(f32_frac) << 29);
+    r64 = call_recip_estimate(f64, 253, fpst);
+    r64_val = float64_val(r64);
+    r64_exp = extract64(r64_val, 52, 11);
+    r64_frac = extract64(r64_val, 0, 52);
+
+    /* result = sign : result_exp<7:0> : fraction<51:29>; */
+    return make_float32(f32_sbit |
+                        (r64_exp & 0xff) << 23 |
+                        extract64(r64_frac, 29, 24));
+}
+
+float64 HELPER(recpe_f64)(float64 input, void *fpstp)
+{
+    float_status *fpst = fpstp;
+    float64 f64 = float64_squash_input_denormal(input, fpst);
+    uint64_t f64_val = float64_val(f64);
+    uint64_t f64_sbit = 0x8000000000000000ULL & f64_val;
+    int64_t f64_exp = extract64(f64_val, 52, 11);
+    float64 r64;
+    uint64_t r64_val;
+    int64_t r64_exp;
+    uint64_t r64_frac;
+
+    /* Deal with any special cases */
+    if (float64_is_any_nan(f64)) {
+        float64 nan = f64;
+        if (float64_is_signaling_nan(f64)) {
+            float_raise(float_flag_invalid, fpst);
+            nan = float64_maybe_silence_nan(f64);
+        }
+        if (fpst->default_nan_mode) {
+            nan =  float64_default_nan;
+        }
+        return nan;
+    } else if (float64_is_infinity(f64)) {
+        return float64_set_sign(float64_zero, float64_is_neg(f64));
+    } else if (float64_is_zero(f64)) {
+        float_raise(float_flag_divbyzero, fpst);
+        return float64_set_sign(float64_infinity, float64_is_neg(f64));
+    } else if ((f64_val & ~(1ULL << 63)) < (1ULL << 50)) {
+        /* Abs(value) < 2.0^-1024 */
+        float_raise(float_flag_overflow | float_flag_inexact, fpst);
+        if (round_to_inf(fpst, f64_sbit)) {
+            return float64_set_sign(float64_infinity, float64_is_neg(f64));
+        } else {
+            return float64_set_sign(float64_maxnorm, float64_is_neg(f64));
+        }
+    } else if (f64_exp >= 1023 && fpst->flush_to_zero) {
+        float_raise(float_flag_underflow, fpst);
+        return float64_set_sign(float64_zero, float64_is_neg(f64));
+    }
 
-    f64 = recip_estimate(f64, env);
+    r64 = call_recip_estimate(f64, 2045, fpst);
+    r64_val = float64_val(r64);
+    r64_exp = extract64(r64_val, 52, 11);
+    r64_frac = extract64(r64_val, 0, 52);
 
-    val32 = sign
-        | ((result_exp & 0xff) << 23)
-        | ((float64_val(f64) >> 29) & 0x7fffff);
-    return make_float32(val32);
+    /* result = sign : result_exp<10:0> : fraction<51:0> */
+    return make_float64(f64_sbit |
+                        ((r64_exp & 0x7ff) << 52) |
+                        r64_frac);
 }
 
 /* The algorithm that must be used to calculate the estimate
@@ -4697,8 +4824,9 @@ float32 HELPER(rsqrte_f32)(float32 a, CPUARMState *env)
     return make_float32(val);
 }
 
-uint32_t HELPER(recpe_u32)(uint32_t a, CPUARMState *env)
+uint32_t HELPER(recpe_u32)(uint32_t a, void *fpstp)
 {
+    float_status *s = fpstp;
     float64 f64;
 
     if ((a & 0x80000000) == 0) {
@@ -4708,7 +4836,7 @@ uint32_t HELPER(recpe_u32)(uint32_t a, CPUARMState *env)
     f64 = make_float64((0x3feULL << 52)
                        | ((int64_t)(a & 0x7fffffff) << 21));
 
-    f64 = recip_estimate (f64, env);
+    f64 = recip_estimate(f64, s);
 
     return 0x80000000 | ((float64_val(f64) >> 21) & 0x7fffffff);
 }