2001-05-15  Andreas Jaeger  <aj@suse.de>

	* sysdeps/ieee754/ldbl-128/s_expm1l.c: New file, contributed by
	Stephen L Moshier <moshier@mediaone.net>.

	* sysdeps/i386/fpu/libm-test-ulps: Adjust for change.

	* math/libm-test.inc: Add comment with ToDo.

	* sysdeps/i386/fpu/e_expl.c: Rewritten to C and using a more
	accurate algorithm.  Patch by Stephen L Moshier <moshier@mediaone.net>.

	* sysdeps/i386/fpu/e_expl.S: Removed.
This commit is contained in:
Andreas Jaeger 2001-05-15 08:00:22 +00:00
parent ef25b29e9a
commit 33996419cd
6 changed files with 109 additions and 54 deletions

View File

@ -1,3 +1,17 @@
2001-05-15 Andreas Jaeger <aj@suse.de>
* sysdeps/ieee754/ldbl-128/s_expm1l.c: New file, contributed by
Stephen L Moshier <moshier@mediaone.net>.
* sysdeps/i386/fpu/libm-test-ulps: Adjust for change.
* math/libm-test.inc: Add comment with ToDo.
* sysdeps/i386/fpu/e_expl.c: Rewritten to C and using a more
accurate algorithm. Patch by Stephen L Moshier <moshier@mediaone.net>.
* sysdeps/i386/fpu/e_expl.S: Removed.
2001-05-14 Roland McGrath <roland@frob.com>
* csu/Makefile ($(objpfx)g$(start-installed-name)): Turn this into a

4
NEWS
View File

@ -6,6 +6,10 @@ Please send GNU C library bug reports using the `glibcbug' script to
<bugs@gnu.org>. Questions and suggestions should be send to
<bug-glibc@gnu.org>.
Version 2.2.4
* Stephen Moshier implemented expm1 for the 128-bit long double format.
Version 2.2.3
* Intel's IA-64 math library is largely integrated. It provides fast and

View File

@ -104,7 +104,12 @@
- Compiler has errors
With e.g. gcc 2.7.2.2 the test for cexp fails because of a compiler error.
*/
To Do: All parameter should be numbers that can be represented as
exact floating point values. Currently some values cannot be represented
exactly and therefore the result is not the expected result.
*/
#ifndef _GNU_SOURCE
# define _GNU_SOURCE

View File

@ -1,43 +0,0 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
*/
#include <machine/asm.h>
RCSID("$NetBSD: $")
/* e^x = 2^(x * log2l(e)) */
ENTRY(__ieee754_expl)
fldt 4(%esp)
/* I added the following ugly construct because expl(+-Inf) resulted
in NaN. The ugliness results from the bright minds at Intel.
For the i686 the code can be written better.
-- drepper@cygnus.com. */
fxam /* Is NaN or +-Inf? */
fstsw %ax
movb $0x45, %dh
andb %ah, %dh
cmpb $0x05, %dh
je 1f /* Is +-Inf, jump. */
fldl2e
fmulp /* x * log2(e) */
fld %st
frndint /* int(x * log2(e)) */
fsubr %st,%st(1) /* fract(x * log2(e)) */
fxch
f2xm1 /* 2^(fract(x * log2(e))) - 1 */
fld1
faddp /* 2^(fract(x * log2(e))) */
fscale /* e^x */
fstp %st(1)
ret
1: testl $0x200, %eax /* Test sign. */
jz 2f /* If positive, jump. */
fstp %st
fldz /* Set result to 0. */
2: ret
END (__ieee754_expl)

75
sysdeps/i386/fpu/e_expl.c Normal file
View File

@ -0,0 +1,75 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
*/
/*
* The 8087 method for the exponential function is to calculate
* exp(x) = 2^(x log2(e))
* after separating integer and fractional parts
* x log2(e) = i + f, |f| <= .5
* 2^i is immediate but f needs to be precise for long double accuracy.
* Suppress range reduction error in computing f by the following.
* Separate x into integer and fractional parts
* x = xi + xf, |xf| <= .5
* Separate log2(e) into the sum of an exact number c0 and small part c1.
* c0 + c1 = log2(e) to extra precision
* Then
* f = (c0 xi - i) + c0 xf + c1 x
* where c0 xi is exact and so also is (c0 xi - i).
* -- moshier@na-net.ornl.gov
*/
static long double c0 = 1.44268798828125L;
static long double c1 = 7.05260771340735992468e-6L;
long double
__ieee754_expl (long double x)
{
long double res, t;
/* I added the following ugly construct because expl(+-Inf) resulted
in NaN. The ugliness results from the bright minds at Intel.
For the i686 the code can be written better.
-- drepper@cygnus.com. */
asm ("fxam\n\t" /* Is NaN or +-Inf? */
"fstsw %%ax\n\t"
"movb $0x45, %%dh\n\t"
"andb %%ah, %%dh\n\t"
"cmpb $0x05, %%dh\n\t"
"je 1f\n\t" /* Is +-Inf, jump. */
"fldl2e\n\t" /* 1 log2(e) */
"fmul %%st(1),%%st\n\t" /* 1 x log2(e) */
"frndint\n\t" /* 1 i */
"fld %%st(1)\n\t" /* 2 x */
"frndint\n\t" /* 2 xi */
"fld %%st(1)\n\t" /* 3 i */
"fldt c0\n\t" /* 4 c0 */
"fld %%st(2)\n\t" /* 5 xi */
"fmul %%st(1),%%st\n\t" /* 5 c0 xi */
"fsubp %%st,%%st(2)\n\t" /* 4 f = c0 xi - i */
"fld %%st(4)\n\t" /* 5 x */
"fsub %%st(3),%%st\n\t" /* 5 xf = x - xi */
"fmulp %%st,%%st(1)\n\t" /* 4 c0 xf */
"faddp %%st,%%st(1)\n\t" /* 3 f = f + c0 xf */
"fldt c1\n\t" /* 4 */
"fmul %%st(4),%%st\n\t" /* 4 c1 * x */
"faddp %%st,%%st(1)\n\t" /* 3 f = f + c1 * x */
"f2xm1\n\t" /* 3 2^(fract(x * log2(e))) - 1 */
"fld1\n\t" /* 4 1.0 */
"faddp\n\t" /* 3 2^(fract(x * log2(e))) */
"fstp %%st(1)\n\t" /* 2 */
"fscale\n\t" /* 2 scale factor is st(1); e^x */
"fstp %%st(1)\n\t" /* 1 */
"fstp %%st(1)\n\t" /* 0 */
"jmp 2f\n\t"
"1:\ttestl $0x200, %%eax\n\t" /* Test sign. */
"jz 2f\n\t" /* If positive, jump. */
"fstp %%st\n\t"
"fldz\n\t" /* Set result to 0. */
"2:\t\n"
: "=t" (res) : "0" (x) : "ax", "dx");
return res;
}

View File

@ -484,13 +484,13 @@ ldouble: 128
# ctan
Test "Real part of: ctan (-2 - 3 i) == 0.0037640256415042482 - 1.0032386273536098014 i":
ildouble: 437
ldouble: 437
ildouble: 439
ldouble: 439
Test "Imaginary part of: ctan (-2 - 3 i) == 0.0037640256415042482 - 1.0032386273536098014 i":
float: 1
ifloat: 1
ildouble: 1
ldouble: 1
ildouble: 2
ldouble: 2
Test "Real part of: ctan (0.7 + 1.2 i) == 0.1720734197630349001 + 0.9544807059989405538 i":
double: 1
float: 1
@ -506,13 +506,13 @@ ldouble: 367
# ctanh
Test "Real part of: ctanh (-2 - 3 i) == -0.9653858790221331242 + 0.0098843750383224937 i":
ildouble: 2
ldouble: 2
ildouble: 5
ldouble: 5
Test "Imaginary part of: ctanh (-2 - 3 i) == -0.9653858790221331242 + 0.0098843750383224937 i":
float: 1
ifloat: 1
ildouble: 23
ldouble: 23
ildouble: 25
ldouble: 25
Test "Real part of: ctanh (0 + pi/4 i) == 0.0 + 1.0 i":
Test "Imaginary part of: ctanh (0 + pi/4 i) == 0.0 + 1.0 i":
float: 1
@ -551,8 +551,8 @@ double: 24
float: 12
idouble: 24
ifloat: 12
ldouble: 4
ildouble: 4
ldouble: 12
ildouble: 12
Test "erfc (9) == 0.41370317465138102381e-36":
ldouble: 36
ildouble: 36