RFC: powerpc: Incorrect results for pow when using FMA

[Adhemerval Zanella <azanella at linux dot vnet dot ibm dot com>]

Hi all,

Checking on a internal issue report, the following testcase is showing
precision errors only for powerpc builds:

--

#include <math.h>
#include <mpfr.h>

int main ()
{
  union Di
  {
    double d;
    unsigned long long i;
  };
  union Di x;
  x.i = 0x553B9BA800000000LL;

  mpfr_t xm;
  mpfr_init2 (xm, 512);
  mpfr_set_d (xm, x.d, MPFR_RNDD);

  mpfr_t squarem;
  mpfr_init2 (squarem, 512);
  mpfr_mul (squarem, xm, xm, MPFR_RNDD);

  double square = x.d * x.d;

  mpfr_printf ("mpfr: x*x   = %.17Re\n", squarem);
  mpfr_printf ("    : x*x   = %.17e\n", square);

  mpfr_t cubem;
  mpfr_init2 (cubem, 512);
  mpfr_mul (cubem, squarem, xm, MPFR_RNDD);

  double cube = square * x.d;

  mpfr_printf ("mpfr: x*x*x = %.17Re\n", cubem);
  mpfr_printf ("    : x*x*x = %.17e\n", cube);

  mpfr_t powm;
  mpfr_init2 (powm, 512);
  mpfr_t onehalfm;
  mpfr_init2 (onehalfm, 512);
  mpfr_set_d (onehalfm, 1.5, MPFR_RNDD);
  mpfr_pow (powm, squarem, onehalfm, MPFR_RNDD);

  double powc = pow (square, 1.5);

  mpfr_printf ("mpfr: pow (x*x, 1.5) = %.17Re\n", powm);
  mpfr_printf ("libc: pow (x*x, 1.5) = %.17e\n", powc);

  return 0;
}

--

With current master I am seeing:

$ ./testrun.sh ./test-mpfr 
mpfr: x*x   = 1.49357829132402765e+205
    : x*x   = 1.49357829132402765e+205
mpfr: x*x*x = 5.77220822056602683e+307
    : x*x*x = 5.77220822056602633e+307
mpfr: pow (x*x, 1.5) = 5.77220822056602683e+307
libc: pow (x*x, 1.5) = 5.77220822056619098e+307

While I would expect the mpfr and libc pow result to have the same 
precision, as for x86_64.

At first I though it was a compiler bug, bug looking into this deeper, 
this is probably a GLIBC build issue.  The difference between the power1() 
-O1 assembler code and -O2 is that the -O2 code makes use of FMA instructions.
This causes a slight difference in values of the second param to the __exp1()
function call.  The __exp1() return value in teh -O1 case returns a invalid
value which forces us to call __slow_pow() leading to the "expected" value. 

In the -O2 case, the __exp1() value is considered valid and we just return 
it as is, leading to the slight difference in values.

If I look at the comments at the top of the source file (sysdeps/ieee754/dbl-64/e_pow.c)
where the power1() function lives, I see this comment:

    /* Assumption: Machine arithmetic operations are performed in              */
    /* round to nearest mode of IEEE 754 standard.                             */

It would seem to me that using FMAs would violate the assumption stated in that comment, 
since some of the internal FAM ops are not rounded before being used.  I do notice that 
if I add the -ffp-contract=off, then we get the "expected" answer.  Is the "correct" 
fix just to get this file compiled with -ffp-contract=off?