Date: Mon, 17 Jun 2002 10:41:17 -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: <3D0E031D.4978E70E@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> <3D0A04A1 DOT 5F275E5C AT cyberoptics DOT com>
 <3D0A53C7 DOT BF363CE0 AT yahoo DOT com>
Reply-To: djgpp-workers AT delorie DOT com

CBFalconer wrote:

> /* =============================================
>  * Check whether (known) integral value is odd */
> bool odd(double val)
> {
>     return ((long)fabs(val) & 1);
> } /* odd */

Unfortunately, this has exactly the same problem as the earlier libc pow(), namely,
that values of "val" greater than (2^31)-1 or less than -(2^31)  or so overflow to
-0x80000000, which tests as even.  Using a cast to "long long" solves the problem
because by the time a "long long" oveflows, all doubles are even ints.  All these
techniques depend on the invalid operation exception being masked, which it is in
DJGPP.  Otherwise an overflow in the cast causes a bomb.  BTW, casting to "long long"
doesn't seem to slow down the routine significantly.

I had earlier chosen to ignore the sign problem for large exponents, because there
even tiny relative errors in the exponent can cause a sign error, which suggests that
in practical computation the results would be garbage anyway.

In fact, there is some justification for returning a NaN if a negative number is
raised to a power larger than 1/DBL_EPSILON, since at that point if the exponent
contains even 1 LSB of rounding error, the sign of the result is indeterminate.
Nonetheless, if one considers the input arguments to pow() to be *exact*, it is not
unreasonable to insist that the sign be computed correctly.

There is a limit to how far this perfectionism can be carried, however, since
eventually one must consider what sign to assign to pow(-1., Inf).  Some people have
argued that infinity is an even integer, which seemed really silly to me.  My libc
pow() function returns a NaN, which seemed like the only responsible thing to do,
since I'm not willing to assume that infinity is an integer at all, let alone an even
integer, and a non-integer power of a negative number is not representable as a
double.

There are still some minor bugs in the libc pow() routine having to do with the signs
of infinite results, but if you take a look at how the various other implementations
of pow() handle these things, you'll see that any caller that depends on these things
being right is treading of very shaky ground.  There are actually about 50 special
cases to resolve in the pow() function, with the treatment of many of them
controversial (such as the example of pow(-1., Inf) above).

-Eric Rudd
rudd AT cyberoptics DOT com