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

[Richard Henderson <rth@twiddle.net>]

On 11/14/2013 06:49 AM, Tom Musta wrote:
> I have reviewed the code and the spec and I cannot see a flaw.  The sequence is
> effectively this:
> 
>   - float64_muladd   - performs proper FMA for 64 bit numbers)
>   - float64_to_float32 - converts to single precision, including proper rounding
>   - float32_to_float64
> 
> The implementation of float64_muladd would seem to provide enough mantissa bits
> for proper handling of the case you describe.  The only rounding occurs in the
> second step.
> 
> I have also done quite a bit of random and targeted random testing using Power
> hardware to produce expected results.  The targeted random tests followed your
> suggestion above: generate AxB + C where abs(exp(A) - exp(B)) = 23 and
> abs(exp(A) - exp(C)) = 46.  Several million test patterns have been generated
> and played back through QEMU without any miscompares in the numerical results.

Here's an example that produces wrong results when rounding to double first.
Replace the portable math.h calls with ppc asm as necessary.


r~


$ cat z.c
#include <stdio.h>
#include <math.h>

float a = 65281;
float b = 257;
float c = 0x1p-29f;

int main()
{
    double dd = fma(a, b, c);
    float d = dd;
    float e = fmaf(a, b, c);

    printf("a = %a\n", a);
    printf("b = %a\n", b);
    printf("c = %a\n", c);
    printf("dd= %a\n", dd);
    printf("d = %a\n", d);
    printf("e = %a\n", e);
    return 0;
}
$ ./a.out
a = 0x1.fe02p+15
b = 0x1.01p+8
c = 0x1p-29
dd= 0x1.000001p+24
d = 0x1p+24
e = 0x1.000002p+24