Date: Fri, 14 Jun 2002 09:58:41 -0500
From: Eric Rudd <rudd AT cyberoptics DOT com>
Subject: Re: [steve AT moshier DOT net: bug in fdlibm e_pow.c]
To: djgpp-workers AT delorie DOT com
Message-id: <3D0A04A1.5F275E5C@cyberoptics.com>
Organization: CyberOptics
MIME-version: 1.0
X-Mailer: Mozilla 4.72 [en] (Win95; U)
Content-type: text/plain; charset=us-ascii
Content-transfer-encoding: 7bit
X-Accept-Language: en,pdf
References: <200206131949 DOT g5DJndk20542 AT envy DOT delorie DOT com>
 <3D09D6D8 DOT 926022C0 AT phekda DOT freeserve DOT co DOT uk>
Reply-To: djgpp-workers AT delorie DOT com

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