1/* ix87 specific implementation of complex exponential function for double.
2   Copyright (C) 1997 Free Software Foundation, Inc.
3   This file is part of the GNU C Library.
4   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
5
6   The GNU C Library is free software; you can redistribute it and/or
7   modify it under the terms of the GNU Lesser General Public
8   License as published by the Free Software Foundation; either
9   version 2.1 of the License, or (at your option) any later version.
10
11   The GNU C Library is distributed in the hope that it will be useful,
12   but WITHOUT ANY WARRANTY; without even the implied warranty of
13   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14   Lesser General Public License for more details.
15
16   You should have received a copy of the GNU Lesser General Public
17   License along with the GNU C Library; if not, write to the Free
18   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
19   02111-1307 USA.  */
20
21#include <sysdep.h>
22
23#ifdef __ELF__
24	.section .rodata
25#else
26	.text
27#endif
28	.align ALIGNARG(4)
29	ASM_TYPE_DIRECTIVE(huge_nan_null_null,@object)
30huge_nan_null_null:
31	.byte 0, 0, 0x80, 0x7f
32	.byte 0, 0, 0xc0, 0x7f
33	.float	0.0
34zero:	.float	0.0
35infinity:
36	.byte 0, 0, 0x80, 0x7f
37	.byte 0, 0, 0xc0, 0x7f
38	.float 0.0
39	.byte 0, 0, 0, 0x80
40	ASM_SIZE_DIRECTIVE(huge_nan_null_null)
41
42	ASM_TYPE_DIRECTIVE(twopi,@object)
43twopi:
44	.byte 0x35, 0xc2, 0x68, 0x21, 0xa2, 0xda, 0xf, 0xc9, 0x1, 0x40
45	.byte 0, 0, 0, 0, 0, 0
46	ASM_SIZE_DIRECTIVE(twopi)
47
48	ASM_TYPE_DIRECTIVE(l2e,@object)
49l2e:
50	.byte 0xbc, 0xf0, 0x17, 0x5c, 0x29, 0x3b, 0xaa, 0xb8, 0xff, 0x3f
51	.byte 0, 0, 0, 0, 0, 0
52	ASM_SIZE_DIRECTIVE(l2e)
53
54	ASM_TYPE_DIRECTIVE(one,@object)
55one:	.double 1.0
56	ASM_SIZE_DIRECTIVE(one)
57
58
59#ifdef PIC
60#define MO(op) op##@GOTOFF(%ecx)
61#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
62#else
63#define MO(op) op
64#define MOX(op,x,f) op(,x,f)
65#endif
66
67	.text
68ENTRY(__cexpf)
69	flds	4(%esp)			/* x */
70	fxam
71	fnstsw
72	flds	8(%esp)			/* y : x */
73#ifdef  PIC
74        call    1f
751:      popl    %ecx
76        addl    $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx
77#endif
78	movb	%ah, %dh
79	andb	$0x45, %ah
80	cmpb	$0x05, %ah
81	je	1f			/* Jump if real part is +-Inf */
82	cmpb	$0x01, %ah
83	je	2f			/* Jump if real part is NaN */
84
85	fxam				/* y : x */
86	fnstsw
87	/* If the imaginary part is not finite we return NaN+i NaN, as
88	   for the case when the real part is NaN.  A test for +-Inf and
89	   NaN would be necessary.  But since we know the stack register
90	   we applied `fxam' to is not empty we can simply use one test.
91	   Check your FPU manual for more information.  */
92	andb	$0x01, %ah
93	cmpb	$0x01, %ah
94	je	20f
95
96	/* We have finite numbers in the real and imaginary part.  Do
97	   the real work now.  */
98	fxch			/* x : y */
99	fldt	MO(l2e)		/* log2(e) : x : y */
100	fmulp			/* x * log2(e) : y */
101	fld	%st		/* x * log2(e) : x * log2(e) : y */
102	frndint			/* int(x * log2(e)) : x * log2(e) : y */
103	fsubr	%st, %st(1)	/* int(x * log2(e)) : frac(x * log2(e)) : y */
104	fxch			/* frac(x * log2(e)) : int(x * log2(e)) : y */
105	f2xm1			/* 2^frac(x * log2(e))-1 : int(x * log2(e)) : y */
106	faddl	MO(one)		/* 2^frac(x * log2(e)) : int(x * log2(e)) : y */
107	fscale			/* e^x : int(x * log2(e)) : y */
108	fst	%st(1)		/* e^x : e^x : y */
109	fxch	%st(2)		/* y : e^x : e^x */
110	fsincos			/* cos(y) : sin(y) : e^x : e^x */
111	fnstsw
112	testl	$0x400, %eax
113	jnz	7f
114	fmulp	%st, %st(3)	/* sin(y) : e^x : e^x * cos(y) */
115	fmulp	%st, %st(1)	/* e^x * sin(y) : e^x * cos(y) */
116	subl	$8, %esp
117	fstps	4(%esp)
118	fstps	(%esp)
119	popl	%eax
120	popl	%edx
121	ret
122
123	/* We have to reduce the argument to fsincos.  */
124	.align ALIGNARG(4)
1257:	fldt	MO(twopi)	/* 2*pi : y : e^x : e^x */
126	fxch			/* y : 2*pi : e^x : e^x */
1278:	fprem1			/* y%(2*pi) : 2*pi : e^x : e^x */
128	fnstsw
129	testl	$0x400, %eax
130	jnz	8b
131	fstp	%st(1)		/* y%(2*pi) : e^x : e^x */
132	fsincos			/* cos(y) : sin(y) : e^x : e^x */
133	fmulp	%st, %st(3)
134	fmulp	%st, %st(1)
135	subl	$8, %esp
136	fstps	4(%esp)
137	fstps	(%esp)
138	popl	%eax
139	popl	%edx
140	ret
141
142	/* The real part is +-inf.  We must make further differences.  */
143	.align ALIGNARG(4)
1441:	fxam			/* y : x */
145	fnstsw
146	movb	%ah, %dl
147	testb	$0x01, %ah	/* See above why 0x01 is usable here.  */
148	jne	3f
149
150
151	/* The real part is +-Inf and the imaginary part is finite.  */
152	andl	$0x245, %edx
153	cmpb	$0x40, %dl	/* Imaginary part == 0?  */
154	je	4f		/* Yes ->  */
155
156	fxch			/* x : y */
157	shrl	$6, %edx
158	fstp	%st(0)		/* y */ /* Drop the real part.  */
159	andl	$8, %edx	/* This puts the sign bit of the real part
160				   in bit 3.  So we can use it to index a
161				   small array to select 0 or Inf.  */
162	fsincos			/* cos(y) : sin(y) */
163	fnstsw
164	testl	$0x0400, %eax
165	jnz	5f
166	fxch
167	ftst
168	fnstsw
169	fstp	%st(0)
170	shll	$23, %eax
171	andl	$0x80000000, %eax
172	orl	MOX(huge_nan_null_null,%edx,1), %eax
173	movl	MOX(huge_nan_null_null,%edx,1), %ecx
174	movl	%eax, %edx
175	ftst
176	fnstsw
177	fstp	%st(0)
178	shll	$23, %eax
179	andl	$0x80000000, %eax
180	orl	%ecx, %eax
181	ret
182	/* We must reduce the argument to fsincos.  */
183	.align ALIGNARG(4)
1845:	fldt	MO(twopi)
185	fxch
1866:	fprem1
187	fnstsw
188	testl	$0x400, %eax
189	jnz	6b
190	fstp	%st(1)
191	fsincos
192	fxch
193	ftst
194	fnstsw
195	fstp	%st(0)
196	shll	$23, %eax
197	andl	$0x80000000, %eax
198	orl	MOX(huge_nan_null_null,%edx,1), %eax
199	movl	MOX(huge_nan_null_null,%edx,1), %ecx
200	movl	%eax, %edx
201	ftst
202	fnstsw
203	fstp	%st(0)
204	shll	$23, %eax
205	andl	$0x80000000, %eax
206	orl	%ecx, %eax
207	ret
208
209	/* The real part is +-Inf and the imaginary part is +-0.  So return
210	   +-Inf+-0i.  */
211	.align ALIGNARG(4)
2124:	subl	$4, %esp
213	fstps	(%esp)
214	shrl	$6, %edx
215	fstp	%st(0)
216	andl	$8, %edx
217	movl	MOX(huge_nan_null_null,%edx,1), %eax
218	popl	%edx
219	ret
220
221	/* The real part is +-Inf, the imaginary is also is not finite.  */
222	.align ALIGNARG(4)
2233:	fstp	%st(0)
224	fstp	%st(0)		/* <empty> */
225	andb	$0x45, %ah
226	andb	$0x47, %dh
227	xorb	%dh, %ah
228	jnz	30f
229	flds	MO(infinity)	/* Raise invalid exception.  */
230	fmuls	MO(zero)
231	fstp	%st(0)
23230:	movl	%edx, %eax
233	shrl	$6, %edx
234	shll	$3, %eax
235	andl	$8, %edx
236	andl	$16, %eax
237	orl	%eax, %edx
238
239	movl	MOX(huge_nan_null_null,%edx,1), %eax
240	movl	MOX(huge_nan_null_null+4,%edx,1), %edx
241	ret
242
243	/* The real part is NaN.  */
244	.align ALIGNARG(4)
24520:	flds	MO(infinity)		/* Raise invalid exception.  */
246	fmuls	MO(zero)
247	fstp	%st(0)
2482:	fstp	%st(0)
249	fstp	%st(0)
250	movl	MO(huge_nan_null_null+4), %eax
251	movl	%eax, %edx
252	ret
253
254END(__cexpf)
255weak_alias (__cexpf, cexpf)
256