On 6 February 2018 at 16:48, Alex Bennée <alex.ben...@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.hender...@linaro.org>
> Signed-off-by: Alex Bennée <alex.ben...@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.mayd...@linaro.org>
thanks
-- PMM
