1/* ix87 specific implementation of exp(x)-1.
2   Copyright (C) 1996, 1997, 2002 Free Software Foundation, Inc.
3   This file is part of the GNU C Library.
4   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
5   Based on code by John C. Bowman <bowman@ipp-garching.mpg.de>.
6   Corrections by H.J. Lu (hjl@gnu.ai.mit.edu), 1997.
7
8   The GNU C Library is free software; you can redistribute it and/or
9   modify it under the terms of the GNU Lesser General Public
10   License as published by the Free Software Foundation; either
11   version 2.1 of the License, or (at your option) any later version.
12
13   The GNU C Library is distributed in the hope that it will be useful,
14   but WITHOUT ANY WARRANTY; without even the implied warranty of
15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
16   Lesser General Public License for more details.
17
18   You should have received a copy of the GNU Lesser General Public
19   License along with the GNU C Library; if not, write to the Free
20   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
21   02111-1307 USA.  */
22
23	/* Using: e^x - 1 = 2^(x * log2(e)) - 1 */
24
25#include <machine/asm.h>
26
27#ifdef __ELF__
28	.section .rodata
29#else
30	.text
31#endif
32	.align ALIGNARG(4)
33	ASM_TYPE_DIRECTIVE(minus1,@object)
34minus1:	.double -1.0
35	ASM_SIZE_DIRECTIVE(minus1)
36	ASM_TYPE_DIRECTIVE(one,@object)
37one:	.double 1.0
38	ASM_SIZE_DIRECTIVE(one)
39	ASM_TYPE_DIRECTIVE(l2e,@object)
40l2e:	.tfloat 1.442695040888963407359924681002
41	ASM_SIZE_DIRECTIVE(l2e)
42
43#ifdef PIC
44#define MO(op) op##@GOTOFF(%edx)
45#else
46#define MO(op) op
47#endif
48
49	.text
50ENTRY(__expm1l)
51	fldt	4(%esp)		// x
52	fxam			// Is NaN or +-Inf?
53	fstsw	%ax
54	movb	$0x45, %ch
55	andb	%ah, %ch
56	cmpb	$0x40, %ch
57	je	3f		// If +-0, jump.
58#ifdef	PIC
59	call	1f
601:	popl	%edx
61	addl	$_GLOBAL_OFFSET_TABLE_+[.-1b], %edx
62#endif
63	cmpb	$0x05, %ch
64	je	2f		// If +-Inf, jump.
65
66	fldt	MO(l2e)		// log2(e) : x
67	fmulp			// log2(e)*x
68	fld	%st		// log2(e)*x : log2(e)*x
69	frndint			// int(log2(e)*x) : log2(e)*x
70	fsubr	%st, %st(1)	// int(log2(e)*x) : fract(log2(e)*x)
71	fxch			// fract(log2(e)*x) : int(log2(e)*x)
72	f2xm1			// 2^fract(log2(e)*x)-1 : int(log2(e)*x)
73	fscale			// 2^(log2(e)*x)-2^int(log2(e)*x) : int(log2(e)*x)
74	fxch			// int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
75	fldl	MO(one)		// 1 : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
76	fscale			// 2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
77	fsubrl	MO(one)		// 1-2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
78	fstp	%st(1)		// 1-2^int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
79	fsubrp	%st, %st(1)	// 2^(log2(e)*x)
80	ret
81
822:	testl	$0x200, %eax	// Test sign.
83	jz	3f		// If positive, jump.
84	fstp	%st
85	fldl	MO(minus1)	// Set result to -1.0.
863:	ret
87END(__expm1l)
88libm_hidden_def (__expm1l)
89weak_alias (__expm1l, expm1l)
90