Bug#182574: libc6: The pow() function gives incorrect results on special cases (nan, inf...)
At Thu, 27 Feb 2003 13:07:40 +0100,
Vincent Lefevre wrote:
> On Thu, Feb 27, 2003 at 12:51:10 +0100, Vincent Lefevre wrote:
> > On Thu, Feb 27, 2003 at 19:59:12 +0900, GOTO Masanori wrote:
> > > Also look at Appendix F.9.4.4. instead of 7.12.7.4-2.
> >
> > OK, I didn't know that they were normalized. I had a quick look at
> >
> > http://www.validlab.com/754R/standards/754.ps
> >
> > (link from <http://www.validlab.com/754R/>), but I didn't see anything
> > concerning transcendental functions. Perhaps an addition of IEC 60559
> > over IEEE 754?
>
> I've just seen:
>
> F.9 Mathematics <math.h>
>
> [#1] This subclause contains specifications of <math.h>
> facilities that are particularly suited for IEC 60559
> implementations.
> [...]
> [#3] Special cases for functions in <math.h> are covered
> directly or indirectly by IEC 60559. The functions that IEC
> 60559 specifies directly are identified in F.3. The other
> functions in <math.h> treat infinities, NaNs, signed zeros,
> subnormals, and (provided the state of the FENV_ACCESS
> pragma is on) the exception flags in a manner consistent
> with the basic arithmetic operations covered by IEC 60559.
>
> So, this seems to be an addition of the ISO C standard.
I'm still investigating this problem, but and I've reached the
conclusion that it seems ISO C99 special rules, so it's not bug.
For example,
http://grouper.ieee.org/groups/754/meeting-materials/2001-07-18-c99.pdf
page 13 contains:
(a) pow(-1, +-inf) = 1
pow(+1, x) = 1 for any x, even a NaN
pow(x, +-0) = 1 for any x, even a NaN
(b) pow(-inf, y) = +0 for any y<0 and not an odd integer
pow(-inf, y) = -inf for y an odd integer > 0
(c) pow(-inf, y) = +inf for y>0 and not an odd integer
So,
Your point is: -> standard says:
1) pow(-inf, nan) = inf -> my test code return it's nan.
2) pow(-inf, 0.5) = inf -> rule (b) can be applied.
pow(-inf, -0.5) = 0 -> rule (c) can be applied.
3) pow(-1, inf) = 1 -> rule (a) can be applied.
pow(-1, -inf) = 1 -> rule (a) can be applied.
Now we've arrived at the question: Is ISO C99 ISO/IEC 60559 floating
point arithmetic really right? There are some discussions about this
issue:
http://www.cs.berkeley.edu/~wkahan/ieee754status/ieee754.ps
Lecture Notes on the Status of IEEE Standard 754 for Binary
Floating-Point Arithmetic by Prof. W. Kahan provides a tour of
some under-appreciated features of IEEE 754 and includes many
examples where they clarify numerical algorithms.
You can see some your interpretations, and you also know much
consideration is needed to change this behavior.
At least all I can say is: it seems that all these special cases is
approved as part of their standard. Glibc math library has a code to
deal with such special cases. It's intentional behavior to match
ISO/IEC 60559.
I think it's already exceeded to discuss on this Debian BTS. If you
think it's bad behavior, please contact to IEEE 754 standard
committee. I would like to close this bug unless you show more
appropriate reason, is it OK?
Regards,
-- gotom
Reply to: