Mail Archives: djgpp-workers/1999/06/08/13:10:41
Eli Zaretskii wrote:
> Perhaps we should discuss the issues raised below on djgpp-workers.
> Eric, if you feel it's a good idea, please redirect the followups as
> appropriate.
>
> On Thu, 3 Jun 1999, DJ Delorie wrote:
>
> > I just got the legal paperwork from Cyberoptics (Eric Rudd's employer)
> > to release some of their libm sources to djgpp. I've put them in
> > ftp://ftp.delorie.com/private/from_erudd.zip
>
> I looked at the code and used these functions with the test programs
> from djtst202.zip. Here's what I found out:
>
> 1) The biggest problem is that asin and acos return NaN for every
> argument. I probably did something stupid to get this,
If it gave you trouble, it'll probably give trouble to dozens of other DJGPP
users, so I had better look into it -- but since they work fine for me, I'm
not sure how to proceed.
> but I
> couldn't find anything that would explain it. Is it possible
> that this is due to some Gas bug (I still have Binutils 2.7
> installed)?
I seem to recall that K. B. Williams was having trouble assembling some of
the routines; in that case it turned out to involve some earlier version of
gas, I think. I'd like to make the routines broadly compatible, but since I
haven't saved those earlier assemblers, I don't have any easy way of
reproducing the problem. I'm using gas v2.8.1.
I have used C-style comments throughout the library, which gives trouble to
earlier versions of gas, unless they are removed by pre-processing by gcc.
The docs have always said one should use C-style comments, so this is
something that is wrong with the earlier versions of gas. Try assembling
the routines with something like
gcc -c asin.s -o asin.o
instead of
as asin.s -o asin.o
In early versions of my code, I had used the # for a comment, since the
earlier assemblers accepted it, but # actually has another meaning in
assembler code, so I got rid of that style of comment.
Let me know if you continue to have trouble.
> 2) Elefunt shows that several functions mis-set errno in some
> cases (see the test programs for details):
>
> - exp2 and exp10 don't set it in both cases which test abnormal
> arguments;
exp2 and exp10 properly set errno when I ran the Elefunt tests. They don't
show up in the .DJL files, however, because they are written to stderr
instead of stdout. Perhaps that is the problem.
> - pow(2, -1.074e3) sets errno to ERANGE whereas Elefunt says it
> shouldn't;
This case produces a denormal. C90 only says in 7.5.1 that "Similarly, a
range error occurs if the result of the function cannot be represented as a
double value." C9x says (in 7.12.1) that "The result underflows if the
magnitude of the mathematical result is so small that the mathematical
result cannot be represented, without extraordinary roundoff error, in an
object of the specified type."
I decided to take the conservative route and flag a range error whenever the
result cannot be represented to full precision (that is, all denormal and
infinite results get flagged with ERANGE). (Actually, pow(2.,-1074.) *can*
be represented exactly, but it is impractical to suppress ERANGE in all such
cases.) I could remove the range error for underflows, but I felt that any
loss of precision due to extreme magnitude of the results ought to be
flagged. Do you think that a different philosophy would be less
problematic?
> - pow(-2, 0) and pow(-2, 2) do NOT set errno, whereas Elefunt
> says they should;
There are some cases, such as pow(0., 0.) or pow(-2., inf), where there is
reasonable debate as to what should be returned, but I can't think of any
reason why we shouldn't have pow(-2.,0.) = 1. and pow(-2., 2.) = 4, with no
error return. I think there's a problem with Elefunt here.
> - sin doesn't set errno at all, while Elefunt expects it for the
> argument of 9.0072e15;
When coding the trig functions, I thought about this case at length, but C90
says "A domain error occurs if an input argument is outside the domain over
which the mathematical function is defined." sin(9.0072E15) is defined in
the mathematical sense, so I decided not to issue a domain error. The
problem with large arguments to sin, cos, and tan is simply a loss of
accuracy, unless the range reduction is done to multiple precision. libm
does this, but a look at their range-reduction code will show that it is no
simple matter. I couldn't think of any practical reason that one should
need such accuracy, since a tiny change in such input arguments makes a
large change in the result. For perfectionists, we'll continue to have libm
(though that library is less accurate in other, more practical cases).
> 3) Paranoia reveals that sqrt(-0.0) returns a NaN,
My personal opinion is that signed zeros are silly, but we've got them. The
rules of IEEE 754 seem to imply that +0 indicates a very small positive
quantity, and -0, a very small negative quantity. Thus, sqrt(-0.) issues a
domain error in my library. I could fix this without any speed penalty.
However, what should sqrt(-0.) return? +0.? -0.? Hmm...
> as does pow(0,0).
> I don't think we want that.
> Also, it might be worth checking
> whether there are other functions that fail for -0.0 but succeed
> for +0.0, and fix that as well.
log(+0.) returns -INF with a range error, but log(-0.) returns NaN with a
domain error, but this seems appropriate, given signed zeros. You may wish
to read the article by John Hauser, "Handling Floating-Point Exceptions in
Numeric Programs:"
http://cch.loria.fr/documentation/IEEE754/ACM/hauser.pdf
In the article, he shows that these decisions are difficult and that no
policy is without its flaws. I'm open to discussion on these matters, but
the choices I've made were generally after much deliberation.
In particular, though the prevailing consensus is that 0^0 should be defined
as 1, my decision to raise EDOM for pow(0., 0.) was based on the fact that
it is mathematically-indeterminate.
> 4) ceil, fabs, floor, and frexp set errno for abnormal arguments,
> but my references indicate that ANSI C89 doesn't require that.
> We are allowed to add error conditions if that is "consistent
> with the mathematical definition of the function", but do we want
> that in these cases?
That's a good question. I don't think that C90 says anything about NaNs,
but C9x makes inconsistent remarks about it: on one hand, in 5.2.4.2.2 it
is stated that
A NaN is an encoding signifying Not-a-Number. A
quiet NaN propagates through almost every arithmetic
operation without raising an exception...
but a NaN is outside the mathematical domain of a function, is it not? It
would actually simplify many routines to simply pass NaNs on through.
> 5) floor, ceil and modf use static variables to store old control
> word, while the old code put them on the stack. This change
> makes functions non-reentrant. Is there a real reason to do
> that?
There are a number of routines that are non-reentrant. If it is important
for them to be reentrant, I will have to fix them.
> 6) hypot will overflow for arguments whose values are larger than
> the square root of LDBL_MAX, even if the result shouldn't
> overflow. The old version didn't have this problem. Is this
> intentional?
hypot is a double routine, not a long double routine. The arguments are
passed as doubles, and it is impossible to have a number as large as
LDBL_MAX in the double format. Thus, I don't see how this can be a
problem. When C9x issues, and support for long double math functions
becomes compulsory, things will become enormously more difficult -- but even
then hypot will not have to change, only hypotl (which hasn't yet been
written).
> 7) Several functions, such as pow and pow10, changed the control
> word in the old version, but now they don't. What was the
> reasons for modifying the control word, and why is it unnecessary
> now?
The question is whether the coprocessor is kept in a "round-to-nearest"
mode. If it is, then the control word does not need to be changed. It is
probably safer to assume nothing, but it is very time-consuming to save the
old control word, reset it, then restore the old control word at exit. Can
we count on the coprocessor being in round-to-nearest mode?
> 8) I understand that we need to keep the old version of modfl, since
> there's no new one, right?
This is a long double (non-ANSI) routine. I don't know why the old libc had
it, since practically no other long double routines were furnished. For
compatibility, we should retain it, however.
> And the same is true for huge_val.c.
Yes.
> 9) We will need to add docs for the following functions:
>
> expm1, cbrt, log1p, exp2, exp10, sincos, log2, powi
>
> The last 5 aren't documented even in libm.info (because libm.a
> doesn't have those functions).
I had started to update the documentation last year, but then things got
bogged down in legal stuff, and meanwhile a new version of DJGPP came out,
with new docs for the math functions. Fortunately, most of my changes were
to indicate ANSI or non-ANSI status, which someone else already did. I had
also intended to indicate the circumstances under which the functions report
domain or range errors, but this is still in flux (see above).
N.B. All but exp10, sincos, powi are specified in the latest draft of C9x.
exp2 is the C9x name for pow2. Though C9x does not define exp10, this
seemed like the appropriate corresponding name for what we now call pow10.
Right now I simply have duplicate labels for them, so they can be called by
either name, and control jumps to the same routine. I had a discussion with
DJ about the necessity for stubs; I have to go back and look that over
again, since I didn't entirely understand the need for them.
I have documentation for sincos and powi at home, which I'll send to you
tomorrow. These functions will probably never make it into standard C, but
I've used them so much myself that I felt they would be useful to others, as
well.
-Eric Rudd
- Raw text -