svn commit: r215237 - head/lib/msun/src

[Bruce Evans]

On Sat, 13 Nov 2010, Alexander Best wrote:

> thank you very much for fixing this long outstanding issue.
> you might want to have a look at [1], where two more issues have been reported.
>
> [1] http://mailman.oakapple.net/pipermail/numeric-interest/2010-September/thread.html

These are either features or not our problem (yet?).  The second issue is
compiler bugs, mostly for VC++.  For gcc, the corresponding bugs are that
gcc is not a C99 compiler, so it doesn't support fenv.  clang has even
worse support for fenv (both have no direct support for fenv, but clang
does more optimizations that are invalid if the environment is not the
default).  fdlibm's problems in this area are that it was written before
fenv existed, so it just assumes that the compiler doesn't do too many
optimizations -- an assumption that is almost but not quite satisifed
with gcc <= 4.2.1 (unless you enable bugs like -ffast-math of course).

The other issue is a (historical?) bug now standardized by C99: from
n1256.pdf:

%          -- pow(-Inf, y) returns +0 for y<0 and not an odd integer.
%          -- pow(-Inf, y) returns +Inf for y>0 and not an odd integer.

So pow(-Inf, 0.5) must return +Inf and pow(-Inf, -0.5) must return +0,
although the bug reporter wants the value to be NaN in both cases,
like it would be for sqrt(-Inf) and 1/sqrt(-Inf), respectively.

pow(-Inf, -0.5) = +0 is reasonable, since it can be interpreted as either
1/sqrt(-Inf) or sqrt(1/-Inf) = sqrt(-0) = 0, and the latter is reasonable
and Standard (sqrt(-0) could be required to be NaN since -0 can be
regarded as being on the branch cut for sqrt(), but that would make -0
too special, especially for real sqrt()).

pow(-Inf, 0.5) = +Inf is not so reasonable, but is consistent with other
(historical?) and now Standard bugs in pow(): essentially, everywhere
that there is a correct special case for the second arg being an odd integer
(so that there is no branch cut involved, and the oddness affects the sign
of the result in the natural way), there is a bogus special case when the
second arg is finite and not an odd integer, and another bogus special case
when the second arge is infinite.  The finite non-odd-integer cases are
generally treated by throwing away all signs so that the result is +Inf
or +0 where it would be +-Inf or +-0 for an odd integer.  Infinite args
are sort of both even and odd, so they should give all of the results
+-Inf or +-0, but this is impossible so the result should be NaN, but
Standards say that it is just +Inf or +0, as results from ignoring the
odd integer singularities.

>From n1156.pdf (not quoting the complementatry odd integer cases):

%          -- pow(+-0, y) returns  +Inf  and  raises  the  divide-by-zero
%             exception for y<0 and not an odd integer.

This is uncontroversial for +0.  However, for -0, it must be interpreted
as 1/pow(-0, -y) and not as 1/pow(1/-0, y) = 1/(-Inf, y) to avoid a
controversial operation on -Inf.

%          -- pow(+-0, y) returns +0 for y>0 and not an odd integer.

Uncontroversial for +0, and less controversial than the above for -0.

%          -- pow(-1, +-Inf) returns 1.

Indirectly related.  This is a new bug in n1156.pdf.  It wasn't in
n869 (not sure if it was in initial-C99), and isn't in fdlibm.  but
this is needed for bug for bug compatibility with pow(-1-eps, +-Inf).
The value of the latter should be NaN, since pow(-1-eps, y) has limits
both +-Inf as y -> +-Inf (strictly the latter should be NaN unless y
is an integer, since a correct complex result is unrepresentable),
but it is +Inf.  But since we fix the sign for the latter, we should
fix the sign for this too.

%          -- pow(x, -Inf) returns +Inf for |x|<1.
%          -- pow(x, -Inf) returns +0 for |x|>1.
%          -- pow(x, +Inf) returns +0 for |x|<1.
%          -- pow(x, +Inf) returns +Inf for |x|>1.

Historical(?) and Standard, but bogus for x < -1 since pow() should only
be defined there when y is an integer and then the result approaches
both +-Inf.  Less bogus for -1 < x <= 0 since then the result (although
it should never be defined :-)) approaches +-0.

%          -- pow(-Inf, y) returns +0 for y<0 and not an odd integer.
%          -- pow(-Inf, y) returns +Inf for y>0 and not an odd integer.

These 2 were quoted earlier.  They are just the rules for +-0 inverted.

cpow() has to deal with pow(x, y) actually existing for negative x, so
it should use projective infinity and not have special cases for +-Inf
or +-0, but presumably it has to be bug for bug compatible with pow()
on the real axis so its special cases can only be more complicated and/or
bogus.

I checked most cases specified in n1156.pdf and found only the following
non-conforming behaviour in FreeBSD:

%C99 rule               -- pow(-Inf, y) returns -0 for y an odd integer < 0.
%fdlibm rule(?)  *	17. -INF ** (anything)  = -0 ** (-anything)

fdlibm error: returns +0 instead of -0.  I'm not sure if I matched the
rules correctly.

%C99 rule               -- pow(-1, +-Inf) returns 1.
%fdlibm rule     *	9.  +-1         ** +-INF is NAN

fdlibm non-error: pow(-1, +-Inf) is NaN, not 1 as specified by C99.
fdlibm non-error: pow(1, +-Inf) is 1 as specified by C99, not NaN as
claimed in the comment.

I know of 1 other similar bug in C99: lgamma(-Inf) should be NaN, but
is Inf.  lgamma() has a very bad singularity at -Inf even when it is
restricted to the real axis.  This bug is missing for tgamma().

Bruce