Fix x86/x86_64 expm1l inaccuracy and exceptions (bugs 13885, 13923).

This commit is contained in:
Joseph Myers 2012-05-07 19:13:08 +00:00
parent 6693d69429
commit 495fd99f3a
9 changed files with 174 additions and 270 deletions

View File

@ -1,3 +1,26 @@
2012-05-07 Joseph Myers <joseph@codesourcery.com>
[BZ #13885]
[BZ #13923]
* sysdeps/i386/fpu/e_expl.S (IEEE754_EXPL): Define conditional on
USE_AS_EXPM1L.
(EXPL_FINITE): Likewise.
(FLDLOG): Likewise.
(IEEE754_EXPL) [USE_AS_EXPM1L]: Support use as expm1l.
* sysdeps/i386/fpu/s_expm1l.S: Define USE_AS_EXPM1L and include
e_expl.S.
* sysdeps/x86_64/fpu/e_expl.S (IEEE754_EXPL): Define conditional on
USE_AS_EXPM1L.
(EXPL_FINITE): Likewise.
(FLDLOG): Likewise.
(IEEE754_EXPL) [USE_AS_EXPM1L]: Support use as expm1l.
* sysdeps/x86_64/fpu/s_expm1l.S: Define USE_AS_EXPM1L and include
e_expl.S.
* math/libm-test.inc (expm1_test): Add more tests. Do not disable
test of -max_value argument for long double.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
2012-05-06 David S. Miller <davem@davemloft.net>
* scripts/data/localplt-sparc-linux-gnu.data: Add '?' markers to

10
NEWS
View File

@ -20,11 +20,11 @@ Version 2.16
13583, 13592, 13618, 13637, 13656, 13658, 13673, 13691, 13695, 13704,
13705, 13706, 13726, 13738, 13739, 13758, 13760, 13761, 13775, 13786,
13787, 13792, 13806, 13824, 13840, 13841, 13844, 13846, 13851, 13852,
13854, 13871, 13872, 13873, 13879, 13883, 13884, 13886, 13892, 13895,
13908, 13910, 13911, 13912, 13913, 13914, 13915, 13916, 13917, 13918,
13919, 13920, 13921, 13922, 13924, 13926, 13927, 13928, 13938, 13941,
13942, 13963, 13967, 13970, 13973, 14027, 14033, 14034, 14040, 14049,
14055, 14064
13854, 13871, 13872, 13873, 13879, 13883, 13884, 13885, 13886, 13892,
13895, 13908, 13910, 13911, 13912, 13913, 13914, 13915, 13916, 13917,
13918, 13919, 13920, 13921, 13922, 13923, 13924, 13926, 13927, 13928,
13938, 13941, 13942, 13963, 13967, 13970, 13973, 14027, 14033, 14034,
14040, 14049, 14055, 14064
* ISO C11 support:

View File

@ -3551,6 +3551,13 @@ expm1_test (void)
TEST_f_f (expm1, 1, M_El - 1.0);
TEST_f_f (expm1, 0.75L, 1.11700001661267466854536981983709561L);
TEST_f_f (expm1, 50.0L, 5.1847055285870724640864533229334853848275e+21L);
#ifndef TEST_FLOAT
TEST_f_f (expm1, 127.0L, 1.4302079958348104463583671072905261080748e+55L);
TEST_f_f (expm1, 500.0L, 1.4035922178528374107397703328409120821806e+217L);
#endif
#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
TEST_f_f (expm1, 11356.25L, 9.05128237311923300051376115753226014206e+4931L);
#endif
@ -3559,9 +3566,7 @@ expm1_test (void)
TEST_f_f (expm1, 100000.0, plus_infty, OVERFLOW_EXCEPTION);
check_int ("errno for expm1(large) == ERANGE", errno, ERANGE, 0, 0, 0);
TEST_f_f (expm1, max_value, plus_infty, OVERFLOW_EXCEPTION);
#ifndef TEST_LDOUBLE /* Bug 13923. */
TEST_f_f (expm1, -max_value, -1);
#endif
END (expm1);
}

View File

@ -28,6 +28,10 @@
# define IEEE754_EXPL __ieee754_exp10l
# define EXPL_FINITE __exp10l_finite
# define FLDLOG fldl2t
#elif defined USE_AS_EXPM1L
# define IEEE754_EXPL __expm1l
# undef EXPL_FINITE
# define FLDLOG fldl2e
#else
# define IEEE754_EXPL __ieee754_expl
# define EXPL_FINITE __expl_finite
@ -69,6 +73,12 @@ csat: .byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x0e, 0x40
.text
ENTRY(IEEE754_EXPL)
#ifdef USE_AS_EXPM1L
movzwl 4+8(%esp), %eax
xorb $0x80, %ah // invert sign bit (now 1 is "positive")
cmpl $0xc006, %eax // is num positive and exp >= 6 (number is >= 128.0)?
jae HIDDEN_JUMPTARGET (__expl) // (if num is denormal, it is at least >= 64.0)
#endif
fldt 4(%esp)
/* I added the following ugly construct because expl(+-Inf) resulted
in NaN. The ugliness results from the bright minds at Intel.
@ -78,7 +88,9 @@ ENTRY(IEEE754_EXPL)
#ifdef PIC
LOAD_PIC_REG (cx)
#endif
#ifndef USE_AS_EXPM1L
movzwl 4+8(%esp), %eax
#endif
andl $0x7fff, %eax
cmpl $0x400d, %eax
jle 3f
@ -96,7 +108,16 @@ ENTRY(IEEE754_EXPL)
andb $2, %ah
jz 3f
fchs
3: FLDLOG /* 1 log2(base) */
3:
#ifdef USE_AS_EXPM1L
/* Test for +-0 as argument. */
fstsw %ax
movb $0x45, %dh
andb %ah, %dh
cmpb $0x40, %dh
je 2f
#endif
FLDLOG /* 1 log2(base) */
fmul %st(1), %st /* 1 x log2(base) */
frndint /* 1 i */
fld %st(1) /* 2 x */
@ -114,17 +135,39 @@ ENTRY(IEEE754_EXPL)
fmul %st(4), %st /* 4 c1 * x */
faddp %st, %st(1) /* 3 f = f + c1 * x */
f2xm1 /* 3 2^(fract(x * log2(base))) - 1 */
#ifdef USE_AS_EXPM1L
fstp %st(1) /* 2 */
fscale /* 2 scale factor is st(1); base^x - 2^i */
fxch /* 2 i */
fld1 /* 3 1.0 */
fscale /* 3 2^i */
fld1 /* 4 1.0 */
fsubrp %st, %st(1) /* 3 2^i - 1.0 */
fstp %st(1) /* 2 */
faddp %st, %st(1) /* 1 base^x - 1.0 */
#else
fld1 /* 4 1.0 */
faddp /* 3 2^(fract(x * log2(base))) */
fstp %st(1) /* 2 */
fscale /* 2 scale factor is st(1); base^x */
fstp %st(1) /* 1 */
#endif
fstp %st(1) /* 0 */
jmp 2f
1: testl $0x200, %eax /* Test sign. */
jz 2f /* If positive, jump. */
fstp %st
#ifdef USE_AS_EXPM1L
fld1
fchs
#else
fldz /* Set result to 0. */
#endif
2: ret
END(IEEE754_EXPL)
#ifdef USE_AS_EXPM1L
libm_hidden_def (__expm1l)
weak_alias (__expm1l, expm1l)
#else
strong_alias (IEEE754_EXPL, EXPL_FINITE)
#endif

View File

@ -1008,10 +1008,6 @@ Test "cos_upward (9) == -0.9111302618846769883682947111811653112463":
ildouble: 1
ldouble: 1
# cosh
Test "cosh (0.75) == 1.29468328467684468784170818539018176":
ildouble: 1
# cosh_downward
Test "cosh_downward (22) == 1792456423.065795780980053377632656584997":
double: 1
@ -1899,29 +1895,20 @@ double: 1
float: 1
idouble: 1
ifloat: 1
ildouble: 4
ldouble: 4
ildouble: 2
ldouble: 2
Test "sinh_downward (23) == 4872401723.124451299966006944252978187305":
double: 1
float: 1
idouble: 1
ifloat: 1
ildouble: 2
ldouble: 2
Test "sinh_downward (24) == 13244561064.92173614705070540368454568168":
float: 1
ifloat: 1
ildouble: 5
ldouble: 5
# sinh_tonearest
Test "sinh_tonearest (22) == 1792456423.065795780701106568345764104225":
ildouble: 3
ldouble: 3
Test "sinh_tonearest (23) == 4872401723.124451299966006944252978187305":
ildouble: 1
ldouble: 1
Test "sinh_tonearest (24) == 13244561064.92173614705070540368454568168":
ildouble: 6
ldouble: 6
ildouble: 2
ldouble: 2
# sinh_towardzero
Test "sinh_towardzero (22) == 1792456423.065795780701106568345764104225":
@ -1929,31 +1916,31 @@ double: 1
float: 1
idouble: 1
ifloat: 1
ildouble: 4
ldouble: 4
ildouble: 2
ldouble: 2
Test "sinh_towardzero (23) == 4872401723.124451299966006944252978187305":
double: 1
float: 1
idouble: 1
ifloat: 1
ildouble: 2
ldouble: 2
Test "sinh_towardzero (24) == 13244561064.92173614705070540368454568168":
float: 1
ifloat: 1
ildouble: 5
ldouble: 5
ildouble: 2
ldouble: 2
# sinh_upward
Test "sinh_upward (22) == 1792456423.065795780701106568345764104225":
ildouble: 16
ldouble: 16
ildouble: 1
ldouble: 1
Test "sinh_upward (23) == 4872401723.124451299966006944252978187305":
ildouble: 27
ldouble: 27
ildouble: 1
ldouble: 1
Test "sinh_upward (24) == 13244561064.92173614705070540368454568168":
double: 1
idouble: 1
ildouble: 7
ldouble: 7
# tan
Test "tan (0x1p16383) == 0.422722393732022337800504160054440141575":
@ -2587,9 +2574,6 @@ ifloat: 1
ildouble: 1
ldouble: 1
Function: "cosh":
ildouble: 1
Function: "cosh_downward":
double: 1
float: 1
@ -2862,26 +2846,22 @@ double: 1
float: 1
idouble: 1
ifloat: 1
ildouble: 5
ldouble: 5
Function: "sinh_tonearest":
ildouble: 6
ldouble: 6
ildouble: 2
ldouble: 2
Function: "sinh_towardzero":
double: 1
float: 1
idouble: 1
ifloat: 1
ildouble: 5
ldouble: 5
ildouble: 2
ldouble: 2
Function: "sinh_upward":
double: 1
idouble: 1
ildouble: 27
ldouble: 27
ildouble: 1
ldouble: 1
Function: "tan":
double: 1

View File

@ -1,89 +1,2 @@
/* ix87 specific implementation of exp(x)-1.
Copyright (C) 1996-1997, 2002, 2005, 2008, 2012 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
Based on code by John C. Bowman <bowman@ipp-garching.mpg.de>.
Corrections by H.J. Lu (hjl@gnu.ai.mit.edu), 1997.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
<http://www.gnu.org/licenses/>. */
/* Using: e^x - 1 = 2^(x * log2(e)) - 1 */
#include <sysdep.h>
#include <machine/asm.h>
.section .rodata
.align ALIGNARG(4)
ASM_TYPE_DIRECTIVE(minus1,@object)
minus1: .double -1.0
ASM_SIZE_DIRECTIVE(minus1)
ASM_TYPE_DIRECTIVE(one,@object)
one: .double 1.0
ASM_SIZE_DIRECTIVE(one)
ASM_TYPE_DIRECTIVE(l2e,@object)
l2e: .tfloat 1.442695040888963407359924681002
ASM_SIZE_DIRECTIVE(l2e)
#ifdef PIC
#define MO(op) op##@GOTOFF(%edx)
#else
#define MO(op) op
#endif
.text
ENTRY(__expm1l)
movzwl 4+8(%esp), %eax // load sign bit and 15-bit exponent
xorb $0x80, %ah // invert sign bit (now 1 is "positive")
cmpl $0xc006, %eax // is num positive and exp >= 6 (number is >= 128.0)?
jae HIDDEN_JUMPTARGET (__expl) // (if num is denormal, it is at least >= 64.0)
fldt 4(%esp) // x
fxam // Is NaN or +-Inf?
fstsw %ax
movb $0x45, %ch
andb %ah, %ch
cmpb $0x40, %ch
je 3f // If +-0, jump.
#ifdef PIC
LOAD_PIC_REG (dx)
#endif
cmpb $0x05, %ch
je 2f // If +-Inf, jump.
fldt MO(l2e) // log2(e) : x
fmulp // log2(e)*x
fld %st // log2(e)*x : log2(e)*x
frndint // int(log2(e)*x) : log2(e)*x
fsubr %st, %st(1) // int(log2(e)*x) : fract(log2(e)*x)
fxch // fract(log2(e)*x) : int(log2(e)*x)
f2xm1 // 2^fract(log2(e)*x)-1 : int(log2(e)*x)
fscale // 2^(log2(e)*x)-2^int(log2(e)*x) : int(log2(e)*x)
fxch // int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fldl MO(one) // 1 : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fscale // 2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fsubrl MO(one) // 1-2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fstp %st(1) // 1-2^int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fsubrp %st, %st(1) // 2^(log2(e)*x)
ret
2: testl $0x200, %eax // Test sign.
jz 3f // If positive, jump.
fstp %st
fldl MO(minus1) // Set result to -1.0.
3: ret
END(__expm1l)
libm_hidden_def (__expm1l)
weak_alias (__expm1l, expm1l)
#define USE_AS_EXPM1L
#include <e_expl.S>

View File

@ -28,6 +28,10 @@
# define IEEE754_EXPL __ieee754_exp10l
# define EXPL_FINITE __exp10l_finite
# define FLDLOG fldl2t
#elif defined USE_AS_EXPM1L
# define IEEE754_EXPL __expm1l
# undef EXPL_FINITE
# define FLDLOG fldl2e
#else
# define IEEE754_EXPL __ieee754_expl
# define EXPL_FINITE __expl_finite
@ -69,13 +73,21 @@ csat: .byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x0e, 0x40
.text
ENTRY(IEEE754_EXPL)
#ifdef USE_AS_EXPM1L
movzwl 8+8(%rsp), %eax
xorb $0x80, %ah // invert sign bit (now 1 is "positive")
cmpl $0xc006, %eax // is num positive and exp >= 6 (number is >= 128.0)?
jae HIDDEN_JUMPTARGET (__expl) // (if num is denormal, it is at least >= 64.0)
#endif
fldt 8(%rsp)
/* 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? */
#ifndef USE_AS_EXPM1L
movzwl 8+8(%rsp), %eax
#endif
andl $0x7fff, %eax
cmpl $0x400d, %eax
jle 3f
@ -93,7 +105,16 @@ ENTRY(IEEE754_EXPL)
andb $2, %ah
jz 3f
fchs
3: FLDLOG /* 1 log2(base) */
3:
#ifdef USE_AS_EXPM1L
/* Test for +-0 as argument. */
fstsw %ax
movb $0x45, %dh
andb %ah, %dh
cmpb $0x40, %dh
je 2f
#endif
FLDLOG /* 1 log2(base) */
fmul %st(1), %st /* 1 x log2(base) */
frndint /* 1 i */
fld %st(1) /* 2 x */
@ -111,17 +132,39 @@ ENTRY(IEEE754_EXPL)
fmul %st(4), %st /* 4 c1 * x */
faddp %st, %st(1) /* 3 f = f + c1 * x */
f2xm1 /* 3 2^(fract(x * log2(base))) - 1 */
#ifdef USE_AS_EXPM1L
fstp %st(1) /* 2 */
fscale /* 2 scale factor is st(1); base^x - 2^i */
fxch /* 2 i */
fld1 /* 3 1.0 */
fscale /* 3 2^i */
fld1 /* 4 1.0 */
fsubrp %st, %st(1) /* 3 2^i - 1.0 */
fstp %st(1) /* 2 */
faddp %st, %st(1) /* 1 base^x - 1.0 */
#else
fld1 /* 4 1.0 */
faddp /* 3 2^(fract(x * log2(base))) */
fstp %st(1) /* 2 */
fscale /* 2 scale factor is st(1); base^x */
fstp %st(1) /* 1 */
#endif
fstp %st(1) /* 0 */
jmp 2f
1: testl $0x200, %eax /* Test sign. */
jz 2f /* If positive, jump. */
fstp %st
#ifdef USE_AS_EXPM1L
fld1
fchs
#else
fldz /* Set result to 0. */
#endif
2: ret
END(IEEE754_EXPL)
#ifdef USE_AS_EXPM1L
libm_hidden_def (__expm1l)
weak_alias (__expm1l, expm1l)
#else
strong_alias (IEEE754_EXPL, EXPL_FINITE)
#endif

View File

@ -1358,6 +1358,9 @@ ifloat: 1
Test "expm1 (11356.25) == 9.05128237311923300051376115753226014206e+4931":
ildouble: 1
ldouble: 1
Test "expm1 (500.0) == 1.4035922178528374107397703328409120821806e+217":
double: 1
idouble: 1
# gamma
Test "gamma (-0.5) == log(2*sqrt(pi))":
@ -1831,62 +1834,47 @@ Test "sincos (pi/6, &sin_res, &cos_res) puts 0.866025403784438646763723170752936
float: 1
ifloat: 1
# sinh
Test "sinh (0x8p-32) == 1.86264514923095703232705808926175479e-9":
ildouble: 1
ldouble: 1
# sinh_downward
Test "sinh_downward (22) == 1792456423.065795780701106568345764104225":
float: 1
ifloat: 1
ildouble: 4
ldouble: 4
ildouble: 2
ldouble: 2
Test "sinh_downward (23) == 4872401723.124451299966006944252978187305":
float: 1
ifloat: 1
ildouble: 2
ldouble: 2
Test "sinh_downward (24) == 13244561064.92173614705070540368454568168":
float: 1
ifloat: 1
ildouble: 5
ldouble: 5
# sinh_tonearest
Test "sinh_tonearest (22) == 1792456423.065795780701106568345764104225":
ildouble: 3
ldouble: 3
Test "sinh_tonearest (23) == 4872401723.124451299966006944252978187305":
ildouble: 1
ldouble: 1
Test "sinh_tonearest (24) == 13244561064.92173614705070540368454568168":
ildouble: 6
ldouble: 6
ildouble: 2
ldouble: 2
# sinh_towardzero
Test "sinh_towardzero (22) == 1792456423.065795780701106568345764104225":
float: 1
ifloat: 1
ildouble: 4
ldouble: 4
ildouble: 2
ldouble: 2
Test "sinh_towardzero (23) == 4872401723.124451299966006944252978187305":
float: 1
ifloat: 1
ildouble: 2
ldouble: 2
Test "sinh_towardzero (24) == 13244561064.92173614705070540368454568168":
float: 1
ifloat: 1
ildouble: 5
ldouble: 5
ildouble: 2
ldouble: 2
# sinh_upward
Test "sinh_upward (22) == 1792456423.065795780701106568345764104225":
ildouble: 16
ldouble: 16
ildouble: 1
ldouble: 1
Test "sinh_upward (23) == 4872401723.124451299966006944252978187305":
ildouble: 27
ldouble: 27
Test "sinh_upward (24) == 13244561064.92173614705070540368454568168":
ildouble: 7
ldouble: 7
ildouble: 1
ldouble: 1
# tan
Test "tan (0x1p16383) == 0.422722393732022337800504160054440141575":
@ -2733,29 +2721,21 @@ ifloat: 1
ildouble: 1
ldouble: 1
Function: "sinh":
ildouble: 1
ldouble: 1
Function: "sinh_downward":
float: 1
ifloat: 1
ildouble: 5
ldouble: 5
Function: "sinh_tonearest":
ildouble: 6
ldouble: 6
ildouble: 2
ldouble: 2
Function: "sinh_towardzero":
float: 1
ifloat: 1
ildouble: 5
ldouble: 5
ildouble: 2
ldouble: 2
Function: "sinh_upward":
ildouble: 27
ldouble: 27
ildouble: 1
ldouble: 1
Function: "tan":
double: 1

View File

@ -1,85 +1,2 @@
/* ix87 specific implementation of exp(x)-1.
Copyright (C) 1996,1997,2001,2002,2008,2009,2012 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
Based on code by John C. Bowman <bowman@ipp-garching.mpg.de>.
Corrections by H.J. Lu (hjl@gnu.ai.mit.edu), 1997.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
<http://www.gnu.org/licenses/>. */
/* Using: e^x - 1 = 2^(x * log2(e)) - 1 */
#include <machine/asm.h>
.section .rodata
.align ALIGNARG(4)
ASM_TYPE_DIRECTIVE(minus1,@object)
minus1: .double -1.0
ASM_SIZE_DIRECTIVE(minus1)
ASM_TYPE_DIRECTIVE(one,@object)
one: .double 1.0
ASM_SIZE_DIRECTIVE(one)
ASM_TYPE_DIRECTIVE(l2e,@object)
l2e: .tfloat 1.442695040888963407359924681002
ASM_SIZE_DIRECTIVE(l2e)
#ifdef PIC
#define MO(op) op##(%rip)
#else
#define MO(op) op
#endif
.text
ENTRY(__expm1l)
movzwl 8+8(%rsp), %eax // load sign bit and 15-bit exponent
xorb $0x80, %ah // invert sign bit (now 1 is "positive")
cmpl $0xc006, %eax // is num positive and exp >= 6 (number is >= 128.0)?
jae __expl // (if num is denormal, it is at least >= 64.0)
fldt 8(%rsp) // x
fxam // Is NaN or +-Inf?
fstsw %ax
movb $0x45, %ch
andb %ah, %ch
cmpb $0x40, %ch
je 3f // If +-0, jump.
cmpb $0x05, %ch
je 2f // If +-Inf, jump.
fldt MO(l2e) // log2(e) : x
fmulp // log2(e)*x
fld %st // log2(e)*x : log2(e)*x
frndint // int(log2(e)*x) : log2(e)*x
fsubr %st, %st(1) // int(log2(e)*x) : fract(log2(e)*x)
fxch // fract(log2(e)*x) : int(log2(e)*x)
f2xm1 // 2^fract(log2(e)*x)-1 : int(log2(e)*x)
fscale // 2^(log2(e)*x)-2^int(log2(e)*x) : int(log2(e)*x)
fxch // int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fldl MO(one) // 1 : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fscale // 2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fsubrl MO(one) // 1-2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fstp %st(1) // 1-2^int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fsubrp %st, %st(1) // 2^(log2(e)*x)-1
ret
2: testl $0x200, %eax // Test sign.
jz 3f // If positive, jump.
fstp %st
fldl MO(minus1) // Set result to -1.0.
3: ret
END(__expm1l)
libm_hidden_def (__expm1l)
weak_alias (__expm1l, expm1l)
#define USE_AS_EXPM1L
#include <e_expl.S>