Message-Id: <200005121101.HAA05022@delorie.com>
From: "Dieter Buerssner" <buers AT gmx DOT de>
To: djgpp-workers AT delorie DOT com
Date: Fri, 12 May 2000 14:09:52 +0200
MIME-Version: 1.0
Content-type: text/plain; charset=US-ASCII
Content-transfer-encoding: 7BIT
Subject: Math functions
X-mailer: Pegasus Mail for Win32 (v3.12b)
Reply-To: djgpp-workers AT delorie DOT com

I have a few comments and questions about the math functions in 
DJGPP. The intention is, that I have contributed a math library to 
the DJGPP project in 1994. It got lost somehow, but I may update it.
It has the long double versions of the most math functions. 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. 

The checking code for the mathfunctions 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.

An example would be the following code:

#include <stdio.h>
#include <math.h>

double my_log(double x)
{
  double y;
  /* Is there a way, to tell gcc, that we need two registers in the 
FPU? This question may be important when the function gets inlined. 
Alternatively, in the input list one could set
: "t" (x), "u" (LN2)
but then, we "loose" one register with fyl2x. Is there a way to tell 
this to gcc?
*/
  __asm__ volatile(
    "fldln2    \n"
    "fldl   %1 \n"
    "fyl2x     \n"
    : "=t" (y)
    : "m" (x));
  return y;
}

int main(void)
{
  double nan;
  volatile short cw = 0x1372; /* exceptions for invalid, div by zero 
and overflow unmasked */
  __asm__ volatile("fldcw %0" : : "m" (cw)); /* Portable to x86 linux 
*/
  nan = log(x); /* No exception, using my_log yields an exception */
  printf("%g\n", nan); /* Exception here, which probably should not 
happen */
  return 0;
}

So, I think, the easier implementation of my_log is preferrable. It 
will always return the same values, but won't set errno.
(my_log should probably be coded as an .S, because of lack of 
documentation for inline assembly with x87-FPU.)

Anyway, is there interest in long double versions of the math 
functions?
Do you prefer, to not raise exceptions?

What I have right now, works with DJGPP 1.x. It is quite fast and 
well tested. To update it for DJGPP 2.x, probably needs only some 
header hacking, to appreciate the new C Standard.

Some functions are a compromise between correctness up to the last 
bit and speed. I.e. the argument reduction code for the trig 
functions is not as sophisticated as the argument reduction code in 
libm, but uses additional precision for intermediate results (from 
Stephen Moshiers cephes), and can be much faster than the 
corresponding code in libm. Similar for powl. I also have functions, 
that are missing right now, like expm1[l] or log1p[l], rint[l] ...

-- 
Regards, Dieter Buerssner