From owner-freebsd-current  Sat Dec 14 03:26:58 1996
Return-Path: <owner-current>
Received: (from root@localhost)
          by freefall.freebsd.org (8.8.4/8.8.4) id DAA22520
          for current-outgoing; Sat, 14 Dec 1996 03:26:58 -0800 (PST)
Received: from godzilla.zeta.org.au (godzilla.zeta.org.au [203.2.228.19])
          by freefall.freebsd.org (8.8.4/8.8.4) with ESMTP id DAA22515
          for <current@FreeBSD.org>; Sat, 14 Dec 1996 03:26:54 -0800 (PST)
Received: (from bde@localhost) by godzilla.zeta.org.au (8.8.3/8.6.9) id WAA02738; Sat, 14 Dec 1996 22:25:06 +1100
Date: Sat, 14 Dec 1996 22:25:06 +1100
From: Bruce Evans <bde@zeta.org.au>
Message-Id: <199612141125.WAA02738@godzilla.zeta.org.au>
To: current@FreeBSD.org, dmaddox@scsn.net
Subject: Re: Something broken in libm
Sender: owner-current@FreeBSD.org
X-Loop: FreeBSD.org
Precedence: bulk

>    There is definitely something broken in libm when HAVE_FPU is
>defined...  More specifically, the problem seems to be in the
>__ieee754_pow() function defined in src/lib/msun/src/e_pow.c.  I
>have not been able to nail down the exact problem, but if anyone
>else is interested --

Try this fix.  __ieee754_pow() only uses two functions that change
when HAVE_FPU is defined, __ieee754_sqrt() and scalbn().
__iee754_sqrt() is to simple to be be broken, but scalbn() left
garbage on the FPU stack.

---
diff -c2 s_scalbn.S~ s_scalbn.S
*** s_scalbn.S~	Fri Aug 19 21:14:30 1994
--- s_scalbn.S	Sat Dec 14 21:36:34 1996
***************
*** 42,44 ****
--- 42,45 ----
  	fldl	4(%esp)
  	fscale
+ 	fstp	%st(1)
  	ret
---

It's surprising that this problem hasn't been noticed before.
scalbn() is also used for log10() and remainders mod pi/2.

Some packages ignore all floating point errors (see fpsetmask(3)).
FPU stack overflow should not be ignored, but unfortunately there
is no way to trap it without also trapping for invalid operands.
When the stack overflows, pushing further operands onto it pushes
NaNs instead of the operands, so nothing much can be computed.
E.g., the following program termintes after 1 iteration instead
of 1e6 iterations:

main()
{
    double d;

    /* Fill up FPU stack. */
    asm("fldz");
    asm("fldz");
    asm("fldz");
    asm("fldz");
    asm("fldz");
    asm("fldz");
    asm("fldz");
    asm("fldz");

    for (d = 0; d < 1e6; )
	++d;
}

Bruce