Re: [Qemu-devel] [V2 PATCH 12/14] target-ppc: VSX Stage 4: Add Scalar SP Fused Multiply-Adds

[Richard Henderson]

On 11/19/2013 11:40 PM, Tom Musta wrote:
> +    /* NOTE: in order to get accurate results, we must first round back */   
>  \
> +    /*       to single precision and use the fused multiply add routine */   
>  \
> +    /*       for 32-bit floats.                                         */   
>  \
> +    float_status tstat = env->fp_status;                                     
>  \
> +    float32 a32 = float64_to_float32(xa.f64[0], &tstat);                     
>  \
> +    float32 b32 = float64_to_float32(b->f64[0], &tstat);                     
>  \
> +    float32 c32 = float64_to_float32(c->f64[0], &tstat);                     
>  \
> +                                                                             
>  \
> +    set_float_exception_flags(0, &tstat);                                    
>  \
> +    float32 t32 = float32_muladd(a32, b32, c32, maddflgs, &tstat);           
>  \

While this will produce correct results for the "normal" use case of correctly
rounded single-precision inputs, the spec says

# Except for xsresp or xsrsqrtesp, any double-precision value can
# be used in single-precision scalar arithmetic operations when
# OE=0 and UE=0.

Thus a more correct implementation would use the full double-precision inputs
while also correctly rounding.  I pointed you at the glibc implementation to
show how that can be done using round-to-zero plus examining the inexact bit.

    float_status tstat = env->fp_status;

    set_float_exception_flags(0, &tstat);
    if (tstat.float_rounding_mode == float_round_nearest_even) {
        /* Avoid double rounding errors by rounding the intermediate
           result to odd.  See
           http://hal.inria.fr/docs/00/08/04/27/PDF/odd-rounding.pdf */
        set_float_rounding_mode(float_round_to_zero, &tstat);
        res = float64_muladd(...);
        res |= (get_float_exception_flags(&tstat) & float_flag_inexact) != 0;
    } else {
        res = float64_muladd(...);
    }
    res = helper_frsp(env, res);

    apply tstat exceptions;


r~