Mail Archives: djgpp-workers/1999/07/04/08:39:16
On Thu, 1 Jul 1999, Eric Rudd wrote:
> No. There are several types of NaN, as you have surmised, but the
> particular NaN reserved by Intel for indefinite results is exactly the bit
> pattern that my functions return (they call it "real indefinite"). Yes,
> the sign bit is set, and I think that's weird, too, but if you check the
> Intel docs, you'll find that the sign bit should indeed be set.
Ugh! This means the value of NaN will be different from what libm
functions return in case of errors, and also from what nan() from libm
returns. Oh well, at least the two bit patterns compare equal, so
users probably won't notice. If nobody else cares, neither do I.
> > It seems to me that if y is finite, x can be Inf and the
> > x87 will produce a finite result, so perhaps we need to modify the
> > checks of x?
>
> You mean, for instance, atan2(1., Inf)? On my platform, this particular
> example returns 0. without error, as it should. I don't see what the
> problem is here. Are we having some more trouble with the assembler?
No problems with assembler here, I get the same results. I added to
the docs the description of the function's behavior for infinite
arguments.
> > - cos suffers a large loss of significant digits near Pi/2: for
> > example cos(-1.571) loses about 40 bits of the mantissa (they come
> > out as zero). The same happens to sin near Pi. Is this some flaw
> > of the x87? If so, perhaps we should instead compute sin(Pi/2 - x)
> > and sin(Pi - x) for such arguments?
>
> The problem is that you need a value of PI that is correct to many more
> than 53 bits in order to produce accurate results under such
> circumstances. The x87 is already doing just such a range reduction, using
> an internal value of PI accurate to 66 bits, but for certain arguments to
> the trig functions (such as the ones you cited), the roundoff errors in the
> range reduction are greatly amplified.
I added a blurb to the docs that is meant for people who might need
the extra accuracy, telling them to link with -lm and pay the price in
slow-down.
> BTW, I noticed that the printf routine prints -0. as +0. (in other
> words, the sign bit is dropped), even when the %+E specifier is used, to
> force a sign to be printed. This doesn't seem right to me. I'd have
> expected %E to print -0. as 0., and %+E to print -0. as -0.
It does print -0.0, but only if the argument is not a real zero
(i.e. some low bits are non-zero, but aren't printed due to format
precision spec). I modified _doprnt so that now it also prints the
minus sign for a negative zero, if the format includes the + flag.
> > - modf doesn't set the sign bit of the integral part if it is zero.
> > For example, -0.81 sets the integral part to +0, not -0. I'm not
> > sure if we should bother.
>
> The positive zero is produced by subtraction of x - fract(x). That's just
> the way it comes out. I'd also lean toward ignoring this little anomaly.
I mentioned this case explcicitly in the docs.
> I'm not sure the present problems with my pow() routine are significant
> enough to hold up the upcoming release, since these behaviors are all ones
> that ANSI does not even specify.
Okay, I agree to postpone the solution of these special cases to the
next release, unless somebody else objects.
- Raw text -