This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Re: [PATCH][ppc] merge in powerpc version of slowpow.c
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: Siddhesh Poyarekar <siddhesh at redhat dot com>
- Cc: <libc-alpha at sourceware dot org>, <munroesj at us dot ibm dot com>, "Ryan S. Arnold"<ryan dot arnold at gmail dot com>
- Date: Thu, 3 Jan 2013 15:44:31 +0000
- Subject: Re: [PATCH][ppc] merge in powerpc version of slowpow.c
- References: <20130103092025.GD5345@spoyarek.pnq.redhat.com>
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