From: rabasse jean-francois <rabasse AT ipno DOT in2p3 DOT fr>
Date: Mon, 9 Sep 1996 16:15:21 +0200 (MET DST)
Message-Id: <199609091415.QAA07831@ipnsub.in2p3.fr>
Content-Type: text
Apparently-To: <djgpp AT delorie DOT com>


Jean-Francois Rabasse, IPN Orsay, France
e-mail rabasse AT ipno DOT in2p3 DOT fr

To:
DJGPP guys & fan-club

Subject: a small computational pbm in djgpp's emu387 package ?

Description:
I build C++ code, with the GCC package, on a PC without math-cop
(80486 SX25), using the djgpp emulator, emu387.dxe. 
A portion of my code computes `acos(x)' values, x being in the range
[+1 ... -1]. I get correct results for acos(+1) to acos(-0.999....)
but acos(-1) returns 0 where PI was expected.

I investigated first the `acos' function source code (in file 
./src/libm/src/e_acos.s). This code computes acos(x) through
atan(sqrt(1 - x^2) / x) using the `fpatan(y / x)' 387 instruction,
and seems Ok.
I looked then at the `fpatan' emulation source code (in file
./src/libemu/emu387.cc) and noticed that the code checks first some
trivial values, nul numerator or denominator of the fpatan(y / x) call.

The reported problem comes from those checks. When calling acos(-1),
the e_acos function computes and calls fpatan(0 / -1) which leaves
the function code through a simple check :
    if (is_zero(st(1)))
    {
        st(1) = CONST_Z;
        ... more code
        return;
    }
The remainder of `fpatan' code does sign checking to setup a `quadrant'
in order to do final value adjustment, and that's why acos(-0.9999..)
returns correct values, near PI.

My comments:
I fixed MY problem by adding a denominator sign check in the above
code :
    if (is_zero(st(1)))
    {
        st(1) = (st(0).sign == SIGN_NEG) ? CONST_PI : CONST_Z;
        ... more code
        return;
    }
but it's not very clever, and doesn't fix ALL problems :

  - the (poor) docs I have, about the FPU instructions tell the fpatan
  expects positive arguments. If so, checking for nul numerator and 
  ignoring the denominator sign is legal and my fix is silly.
  (By the way, a similar check is done for nul denominator to return
  PI/2 value whatever the numerator sign is, so fpatan(+x / 0) will
  return PI/2 but fpatan(-x / 0) will return the same, not -PI/2 !)

  - implementing a fpatan function accepting only positive arguments
  may break the acos function implementation. From a math point of
  view, the atan() function is supposed to return a value in the
  range [-PI/2 ... +PI/2], and the acos function is supposed to return
  a value in the range [0 ... PI].
  An `acos' implementation should be atan(sqrt(1 - x^2) / x) for
  non negative x, and (atan(sqrt(1 - x^2) / x) + PI) for negative x,
  provided the `atan' implementation returns values in range 
  [-PI/2 ... +PI/2] and accepts negative arguments.

Conclusion:
I ignore if this problem is known and have yet been fixed.
If yes, please burn this mail.

If no, it's probably not worth a fix. All CPU's after the 486 SX have
now a FPU chip, and warning potential users of the emulib, with old
machines like mine, may be suffisant.

Sincerely, JF-R