Date: Fri, 12 May 2000 11:55:36 -0500
From: Eric Rudd <rudd AT cyberoptics DOT com>
Subject: Re: Math functions
To: djgpp-workers AT delorie DOT com
Message-id: <391C3788.70372F70@cyberoptics.com>
Organization: CyberOptics
MIME-version: 1.0
X-Mailer: Mozilla 4.72 [en] (Win95; U)
Content-type: text/plain; charset=us-ascii
Content-transfer-encoding: 7bit
X-Accept-Language: en
References: <200005121101 DOT HAA05022 AT delorie DOT com>
Reply-To: djgpp-workers AT delorie DOT com

Dieter Buerssner wrote:

> While looking at the current math functions in libc, I see, that it
> tries hard to find special arguments and domain errors, and sets errno.
> The setting of errno seems not to be required by the latest C Standard.
> Is the setting of errno really wanted for DJGPP? To me, it seems
> redundant, because the error can already be detected by the return
> value.

It's not as simple as that.  An implementation is not required to set
errno under C99, and initially I didn't want to set it, either, but Eli
convinced me that setting it under error conditions would be a good
idea.  It turned out not to slow down the routines much, and helps in
cases like this:

   for (i=0; i<n; i++) {
      a[i] = log(b[i]);
   }
   if (errno != 0) ...

where one can save the time of checking each individual result.
Furthermore, an implementation is not required to return NaN if there's
an error, so in the absence of proper errno support, one has to put in
more complicated tests, which should really have gone inside the
function.  For instance, log(0.) returns -Inf instead, but it's a range
error in my implementation.

There is some defense for the point of view that a portable program
shouldn't count on implementation-defined features, and that therefore
the implementation itself shouldn't include any features not *required*
by the standard.  Unfortunately, the C standards don't say much about
accuracy, either, so must a portable program not even count on accurate
functions?  It still boils down to a quality-of-implementation issue, and
errno contributes to a high-quality implementation.

> The checking code for the math functions silently loads a NaN into the
> top of the FPU-stack, when a domain error occures. This has IHMO a
> significant disadvantage. When setting up the FPU control word, to get
> exceptions for (say) overflow, div by zero and invalid (which I find
> useful while debugging), the functions won't raise that exception. So
> the error may be detected much later, than the point were it occured,
> or even slip by undetected.

This behavior is required by C90:

     7.5.1  Treatment of error conditions

        The behavior of each of these functions is defined for all
     representable values of its input arguments.  Each function
     shall execute as if it were a single operation, without
     generating any externally visible exceptions.

> Anyway, is there interest in long double versions of the math
> functions?

I think they'll be needed for C99.  The interesting question is whether
the coprocessor accuracy should be considered adequate or not.  Long
doubles have a 64-bit significand, which is the same as the internal
format of the x87, but it's not easy to deliver correctly-rounded results
at that precision for more than a few of the functions (such as sqrt).
The range reduction for exp(), sinh(), cosh(), and the trig functions
will be particularly difficult.

If a program uses long doubles, it is presumably very concerned about
accuracy, so I would lean toward slower, uncompromisingly-accurate
implementations of the long double functions.  I'd want to make sure that
such implementations did not get linked into programs not using long
doubles, however.

> Do you prefer, to not raise exceptions?

See above.

> I also have functions, that are missing right now, like expm1[l] or
> log1p[l], rint[l] ...

libc already has expm1() and log1p(), but not rint.

I contributed the "double" math functions in libc, and am willing to
update them to C99, if and when the standard is ratified here in the US,
but I'm less interested in the long double functions.  If you would like
to work on long double routines, that's fine with me.

One issue for me is that C99 requires certain behaviors that I consider
mathematically indefensible, such as atan2(0.,0.) returning 0. without
any indication of error.  I would hate to update my routines to behave in
that way, even though it would make them faster and simpler.  Certain
people (myself included) submitted formal public comments about this, and
I understand that one of the other people has pursued it a little more
aggressively, which is delaying ratification here in the US.  I don't
know what the outcome will be, but I agree with his objections.  I'd like
to wait until things settle down before proceeding to update the
implementations.

There are a number of higher mathematical functions that will need to be
implemented for C99, such as lgamma, tgamma, erf, erfc, etc.  I don't
know if we want to adopt the routines in libm or if some improvement can
be made in them.  I have a rational minimax program that would be useful
in cranking out coefficients for these functions, but it would need to be
upgraded to use a multiple-precision math package in order to be useful
for double and long double routines.  There are also complex routines
that will be needed, and there are, as I recall, some non-trivial issues
with maintaining accuracy there.

-Eric Rudd
rudd AT cyberoptics DOT com