Mail Archives: djgpp-workers/2002/06/14/12:30:19
Richard Dawe wrote:
> After applying the patch I get:
>
> bash-2.04$ ./pow-bug
> pow(1.0000000000000002e+00 4.5035996273704970e+15) = 2.7182818284590455e+00
> pow(-1.0000000000000002e+00 4.5035996273704970e+15) = 2.7182818284590455e+00
> pow(9.9999999999999978e-01 4.5035996273704970e+15) = 3.6787944117144222e-01
> pow(-9.9999999999999978e-01 4.5035996273704970e+15) = 3.6787944117144222e-01
> bash-2.04$ ./pow-bug-m
> pow(1.0000000000000002e+00 4.5035996273704970e+15) = 2.7182818284590455e+00
> pow(-1.0000000000000002e+00 4.5035996273704970e+15) = -2.7182818284590455e+00
> pow(9.9999999999999978e-01 4.5035996273704970e+15) = 3.6787944117144222e-01
> pow(-9.9999999999999978e-01 4.5035996273704970e+15) = -3.6787944117144222e-01
>
> Note that the results are not consistent - the sign differs in libc vs. libm.
My libc routine is wrong -- see below.
> Now, if I can still remember my maths correctly 8) :
>
> a ** b == exp(b * ln(a))
>
> So if a is negative then ln(a) gives a complex number. How does a function
> that deals with real numbers cope with the imaginary part of that?
That formula is mathematically correct, but impractical to use for negative a.
In practice, for negative a one computes
a^b = (sign(a)*|a|)^b = (sign(a))^b * exp(b*ln(|a|))
The exp() term is as easy to compute when a is negative as when it's positive.
The first term simply corrects the sign: for a>=0, the result is never negative;
for a<0, if b is even, the result is positive; if b is odd, the result is
negative; if b is non-integer, pow() returns NaN, since the result is complex.
Since I wrote the pow() function in libc, I figured I was in as good a position
as anyone to figure out the cause of the discrepancy. I took a look at my code,
and very quickly saw the reason: I had used the "fistl" instruction to write b
out to memory as a long int. I then tested its evenness by and'ing it with 1.
This test works as long as the number is within the range of a long int, but
fails for really big numbers such as Moshier was using, since fistl clips (or
saturates) out-of-range numbers, instead of writing them out modulo 2^32, as I
wish it would. I'll take a look at my code, and see if there's a reasonably
efficient test of evenness for big integers. There's a "fistq" instruction,
which writes out 64-bit integers and would presumably solve the problem, but it
will be slower. The reason 64-bit integers should solve the problem is that any
doubles > 2^64 in absolute value are always even integers, since their
significands are only 53 bits.
-Eric Rudd
rudd AT cyberoptics DOT com
- Raw text -