This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: [PATCH][ppc] merge in powerpc version of slowpow.c


On Thu, 3 Jan 2013, Siddhesh Poyarekar wrote:

> +#ifndef NO_LONG_DOUBLE
> +  long double ldw, ldz, ldpp;
> +  static const long double ldeps = 0x4.0p-96;
> +
> +  /* Compute pow as long double, 106 bits */
> +  ldz = __ieee754_logl ((long double) x);
> +  ldw = (long double) y *ldz;
> +  ldpp = __ieee754_expl (ldw);

The logic for 0x4.0p-96 in the IBM long double case would be that up to 10 
bits of the 106-bit value of ldw might end up as exponent bits rather than 
manitssa bits in the value of ldpp, so that the mantissa of ldpp is only 
accurate to 96 bits.  There are two problems there even for IBM long 
double:

* Justifying the value "4" as big enough to cover all possible error 
accumulation in the calculation requires a detailed error analysis of 
logl, multiplication, expl and the addition / subtraction used in testing 
whether accumulated errors could have affected the rounding (recall that 
multiplication and addition are *not* 0.5ulp accurate for IBM long 
double); it's quite possible the errors could accumulate to more than 4 
parts in 2**96 for some inputs.

* The above argument is for a *relative* error of 0x4.0p-96 (or some value 
of that order of magnitude), but the logic in the code is treating it as 
an absolute error (meaning that it will wrongly deduce for some large 
results that they must be correctly rounded from the long double 
computation, when they may not be).

So I don't think the IBM long double code being used on Power is correct 
as-is, because it appears to be using 0x4.0p-96 as an absolute error 
value; computing e.g. ldpp * (1.0L + ldeps) would seem more appropriate.  
And for x86 long double with 64-bit mantissa, this approproach isn't 
useful at all, since you can't get more than about 54 bits precision from 
a calculation using 64-bit mantissa internally, which means you never have 
enough information to know what the correctly rounded value was.

-- 
Joseph S. Myers
joseph@codesourcery.com


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]