Mail Archives: djgpp-workers/2002/06/17/12:00:45
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
- Raw text -