Re: [Qemu-devel] [PATCH v3 22/31] arm/helper.c: re-factor recpe and add recepe_f16

[Richard Henderson]

On 02/23/2018 07:36 AM, Alex Bennée wrote:
> It looks like the ARM ARM has simplified the pseudo code for the
> calculation which is done on a fixed point 9 bit integer maths. So
> while adding f16 we can also clean this up to be a little less heavy
> on the floating point and just return the fractional part and leave
> the calle's to do the final packing of the result.
> 
> Signed-off-by: Alex Bennée 
> 
> ---
> v3
>   - fix comment 2.0^-16
>   - f16_exp >= 29 (biased)
>   - remove confusing comment about fpst
> ---

Reviewed-by: Richard Henderson 


r~

[Qemu-devel] [PATCH v3 22/31] arm/helper.c: re-factor recpe and add recepe_f16

[Alex Bennée]

It looks like the ARM ARM has simplified the pseudo code for the
calculation which is done on a fixed point 9 bit integer maths. So
while adding f16 we can also clean this up to be a little less heavy
on the floating point and just return the fractional part and leave
the calle's to do the final packing of the result.

Signed-off-by: Alex Bennée 

---
v3
  - fix comment 2.0^-16
  - f16_exp >= 29 (biased)
  - remove confusing comment about fpst
---
 target/arm/helper.c | 224 +---
 target/arm/helper.h |   1 +
 2 files changed, 128 insertions(+), 97 deletions(-)

diff --git a/target/arm/helper.c b/target/arm/helper.c
index 6e3dadb754..e2d0ff0b4c 100644
--- a/target/arm/helper.c
+++ b/target/arm/helper.c
@@ -11523,80 +11523,75 @@ float32 HELPER(rsqrts_f32)(float32 a, float32 b, 
CPUARMState *env)
  * int->float conversions at run-time.  */
 #define float64_256 make_float64(0x4070LL)
 #define float64_512 make_float64(0x4080LL)
+#define float16_maxnorm make_float16(0x7bff)
 #define float32_maxnorm make_float32(0x7f7f)
 #define float64_maxnorm make_float64(0x7fefLL)
 
 /* Reciprocal functions
  *
  * The algorithm that must be used to calculate the estimate
- * is specified by the ARM ARM, see FPRecipEstimate()
+ * is specified by the ARM ARM, see FPRecipEstimate()/RecipEstimate
  */
 
-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 = *real_fp_status;
-float_status *s = _status;
-/* q = (int)(a * 512.0) */
-float64 q = float64_mul(float64_512, a, s);
-int64_t q_int = float64_to_int64_round_to_zero(q, s);
-
-/* r = 1.0 / (((double)q + 0.5) / 512.0) */
-q = int64_to_float64(q_int, s);
-q = float64_add(q, float64_half, s);
-q = float64_div(q, float64_512, s);
-q = float64_div(float64_one, q, s);
-
-/* s = (int)(256.0 * r + 0.5) */
-q = float64_mul(q, float64_256, s);
-q = float64_add(q, float64_half, s);
-q_int = float64_to_int64_round_to_zero(q, s);
+/* See RecipEstimate()
+ *
+ * input is a 9 bit fixed point number
+ * input range 256 .. 511 for a number from 0.5 <= x < 1.0.
+ * result range 256 .. 511 for a number from 1.0 to 511/256.
+ */
 
-/* return (double)s / 256.0 */
-return float64_div(int64_to_float64(q_int, s), float64_256, s);
+static int recip_estimate(int input)
+{
+int a, b, r;
+assert(256 <= input && input < 512);
+a = (input * 2) + 1;
+b = (1 << 19) / a;
+r = (b + 1) >> 1;
+assert(256 <= r && r < 512);
+return r;
 }
 
-/* Common wrapper to call recip_estimate */
-static float64 call_recip_estimate(float64 num, int off, float_status *fpst)
+/*
+ * Common wrapper to call recip_estimate
+ *
+ * The parameters are exponent and 64 bit fraction (without implicit
+ * bit) where the binary point is nominally at bit 52. Returns a
+ * float64 which can then be rounded to the appropriate size by the
+ * callee.
+ */
+
+static uint64_t call_recip_estimate(int *exp, int exp_off, uint64_t frac)
 {
-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;
+uint32_t scaled, estimate;
+uint64_t result_frac;
+int result_exp;
 
-/* Generate the scaled number for the estimate function */
-if (exp == 0) {
+/* Handle sub-normals */
+if (*exp == 0) {
 if (extract64(frac, 51, 1) == 0) {
-exp = -1;
-frac = extract64(frac, 0, 50) << 2;
+*exp = -1;
+frac <<= 2;
 } else {
-frac = extract64(frac, 0, 51) << 1;
+frac <<= 1;
 }
 }
 
-/* scaled = '0' : '010' : fraction<51:44> : Zeros(44); */
-scaled = make_float64((0x3feULL << 52)
-  | extract64(frac, 44, 8) << 44);
-
-estimate = recip_estimate(scaled, fpst);
+/* scaled = UInt('1':fraction<51:44>) */
+scaled = deposit32(1 << 8, 0, 8, extract64(frac, 44, 8));
+estimate = recip_estimate(scaled);
 
-/* Build new result */
-val64 = float64_val(estimate);
-sbit = 0x8000ULL & 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;
+result_exp = exp_off - *exp;
+result_frac = deposit64(0, 44, 8, estimate);
+if (result_exp == 0) {
+result_frac = deposit64(result_frac >> 1, 51, 1, 1);
+} else if (result_exp == -1) {
+result_frac = deposit64(result_frac >> 2, 50, 2, 1);
+result_exp = 0;
 }
 
-return make_float64(sbit | (exp << 52) | frac);
+*exp = result_exp;