Message-Id: <200005131046.GAA30186@delorie.com>
From: "Dieter Buerssner" <buers AT gmx DOT de>
To: djgpp-workers AT delorie DOT com
Date: Sat, 13 May 2000 13:54:37 +0200
MIME-Version: 1.0
Content-type: text/plain; charset=US-ASCII
Content-transfer-encoding: 7BIT
Subject: Re: Math functions
In-reply-to: <391C7D31.838DE634@cyberoptics.com>
X-mailer: Pegasus Mail for Win32 (v3.12b)
Reply-To: djgpp-workers AT delorie DOT com
Errors-To: nobody AT delorie DOT com
X-Mailing-List: djgpp-workers AT delorie DOT com
X-Unsubscribes-To: listserv AT delorie DOT com
Precedence: bulk

On 12 May 00, at 16:52, Eric Rudd wrote:

> Dieter Buerssner wrote:
> > I think the range reduction for expl() is not that difficult, because of the
> > limited range of the argument. When expl() and expm1l() are implemented,
> > sinhl(), coshl() and tanhl() are not that difficult either, when an maximum
> > errror of say 2ULPs = 2^-62 is acceptable.
> 
> I haven't looked at Moshier's code lately, but the inherent problem with range
> reduction for functions like expl() is that the range reduction needs to be
> done to higher precision than that of the final result, because of the
> magnfication of the relative error of the argument.  It turns out that for
> expl(), with an output significand of 64 bits, and an exponent range of  +/-
> 16384 (or 15 bits), you need about 15 bits of extra precision in the range
> reduction, or about 64+15=79.  The 64-bit precision of the coprocessor was
> adequate for range reduction in a double routine, but not long double.

I am using at least 16 additional bits. The following code was 
inspired by C code of Moshier, from which I learned the extended 
precision tricks. It needs just 3 fmuls, 2 fsubs and a few flds more, 
than the code without additional precision.

/* expl for arguments for which abs(log2(e)*x) < 2^48
   Note, that the result will overflow or underflow long before that.
   The argument reduction is done with more than 16 additional bits,
   which seems enough over the whole range. (The commented constants
   are for 32 additional bits, which means, that abs(log2(2) * x)
   must be smaller 2^32. The algorithm should not depend on the
   rounding mode, but only round to nearest was tested.

      e^x = 2^(x*log2(e))
          = 2^(Integer(x*log2(e)) * 2^(x*log2(e) - Integer(x*log2(e)))
      x*log2(e) - Integer(x*log2(e))
          = (x - Integer(x*log2(e)) * log(2)) * log2(e)

   The subtraction in the last formula is done with additional precision.
   LN2L_HIGH * Integer(x*log2(e)) will always yield an exact result (no
   rounding involved).

   The constants are coded in hex, to not depend on correct
   dec -> bin conversion of the compiler/assembler.

   -- Dieter Buerssner
*/

	.data
/* LOG2L_HIGH + LOG2L_LOW = 6.931471805599453094172321214581765680755E-1 */
LN2L_HIGH:
	.long 0x0,0xb1720000,0x3ffe /* only 16 high mantissa bits set */
LN2L_LOW:
	.long 0xcd5e4f1e,0xbfbe8e7b,0x3feb /* what remains */

/* 32 excess bit constants
 0x0,0xb17217f8,0x3ffe
 0xd87131a0,0xb8c21950,0xbfdc */

	.text
	.globl ___expl
___expl:
	pushl %ebp
	movl %esp,%ebp
	fldt 8(%ebp)
	fldl2e
	fmul %st(1),%st
	frndint             /* ix=Integer(log2(e)*x), x */
	fldt LN2L_HIGH
	fmul %st(1),%st
	fsubrp %st,%st(2)   /* x -= ix * LN2L_HIGH */
	fldt LN2L_LOW
	fmul %st(1),%st
	fsubrp %st,%st(2)   /* x -= ix * LN2L_LOW */
	fldl2e
	fmulp %st,%st(2)    /* x *= log2(e) */
	fxch %st(1)
	f2xm1               /* 2^x - 1 */
	fld1
	faddp %st,%st(1)    /* 2^x */
	fscale              /* 2^ix * 2^x */
	fstp %st(1)
	leave
	ret

Supplemented by some C code:

#include <errno.h>

#define maxlogl  1.1356523406294143949492E4L
#define minlogl -1.1355137111933024058873E4L

static float one = 1.0, zero=0.0;
static long double tiny = 1e-3000L;

extern long double __expl(long double);

long double expl(long double x)
{
  /* Could be changed to some checks of x at the bit level.
     The exact values of maxlogl and minlogl are not that important.
     __expl will overflow or underflow anyway, but errno would
     be lost. */
  if (x != x)
    return x+x; /* NaN */
  if(x > maxlogl)
  {
    errno = ERANGE;
    return one/zero; /* trigger overflow */
  }
  if(x < minlogl)
    return tiny*tiny; /* trigger underflow */
  return __expl(x);
}
 
> > Sin(), cos() and tan() are more difficult.

> It strikes me that it's difficult for any argument where fabs(x) >> 1.   K. B.
> Williams pointed me to some technical articles showing how to compute a
> correctly-rounded double result for any double argument:
> 
>     http://www.validgh.com/arg.ps
> 
> but it seemed like an unreasonable amount of extra execution time for
> "workhorse" code that might well be used where speed is critical. 

I totally agree. (I have read that article quite some time ago.) Also 
it's beyond my skills, to write the argument reduction code, with 
2^14 additional bits of precision. Also, when anybody calls 
sinl(1e1000L), he is probably already in trouble. The closest 
representable values to 1e1000L will be many periods of sin apart. I 
really cannot imagine any application, where this could be useful.

> > The most difficult is IHMO powl(), and I don't have a good solution
> > for that.
> 
> powl() is definitely one of the most difficult functions to compute.  If one
> carries internal results to 80 bits, normal arguments can be handled without
> undue difficulty, 

You mean 80 bits in the mantissa, yes?
Do you have code for 80 bit log and exp?

> but there are about 50 exception conditions that must be
> handled -- powl(inf, 0.), powl(0., NaN), powl(-1., inf), and so forth.  (Let's
> see... is infinity an integer? ;-) 

The list in my C99 draft is horribly long ...

-- 
Regards, Dieter Buerssner