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, 0, 0, 0, 0, 0xf0, 0x7f
32	.byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
33	.double	0.0
34zero:	.double	0.0
35infinity:
36	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
37	.byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
38	.double 0.0
39	.byte 0, 0, 0, 0, 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(__cexpl)
69	fldt	8(%esp)			/* x */
70	fxam
71	fnstsw
72	fldt	20(%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	movl	4(%esp), %eax		/* Pointer to memory for result.  */
117	fstpt	12(%eax)
118	fstpt	(%eax)
119	ret	$4
120
121	/* We have to reduce the argument to fsincos.  */
122	.align ALIGNARG(4)
1237:	fldt	MO(twopi)	/* 2*pi : y : e^x : e^x */
124	fxch			/* y : 2*pi : e^x : e^x */
1258:	fprem1			/* y%(2*pi) : 2*pi : e^x : e^x */
126	fnstsw
127	testl	$0x400, %eax
128	jnz	8b
129	fstp	%st(1)		/* y%(2*pi) : e^x : e^x */
130	fsincos			/* cos(y) : sin(y) : e^x : e^x */
131	fmulp	%st, %st(3)
132	fmulp	%st, %st(1)
133	movl	4(%esp), %eax		/* Pointer to memory for result.  */
134	fstpt	12(%eax)
135	fstpt	(%eax)
136	ret	$4
137
138	/* The real part is +-inf.  We must make further differences.  */
139	.align ALIGNARG(4)
1401:	fxam			/* y : x */
141	fnstsw
142	movb	%ah, %dl
143	testb	$0x01, %ah	/* See above why 0x01 is usable here.  */
144	jne	3f
145
146
147	/* The real part is +-Inf and the imaginary part is finite.  */
148	andl	$0x245, %edx
149	cmpb	$0x40, %dl	/* Imaginary part == 0?  */
150	je	4f		/* Yes ->  */
151
152	fxch			/* x : y */
153	shrl	$5, %edx
154	fstp	%st(0)		/* y */ /* Drop the real part.  */
155	andl	$16, %edx	/* This puts the sign bit of the real part
156				   in bit 4.  So we can use it to index a
157				   small array to select 0 or Inf.  */
158	fsincos			/* cos(y) : sin(y) */
159	fnstsw
160	testl	$0x0400, %eax
161	jnz	5f
162	fldl	MOX(huge_nan_null_null,%edx,1)
163	movl	4(%esp), %edx		/* Pointer to memory for result.  */
164	fld	%st
165	fstpt	12(%edx)
166	fstpt	(%edx)
167	ftst
168	fnstsw
169	shll	$7, %eax
170	andl	$0x8000, %eax
171	orl	%eax, 8(%edx)
172	fstp	%st(0)
173	ftst
174	fnstsw
175	shll	$7, %eax
176	andl	$0x8000, %eax
177	orl	%eax, 20(%edx)
178	fstp	%st(0)
179	ret	$4
180	/* We must reduce the argument to fsincos.  */
181	.align ALIGNARG(4)
1825:	fldt	MO(twopi)
183	fxch
1846:	fprem1
185	fnstsw
186	testl	$0x400, %eax
187	jnz	6b
188	fstp	%st(1)
189	fsincos
190	fldl	MOX(huge_nan_null_null,%edx,1)
191	movl	4(%esp), %edx		/* Pointer to memory for result.  */
192	fld	%st
193	fstpt	12(%edx)
194	fstpt	(%edx)
195	ftst
196	fnstsw
197	shll	$7, %eax
198	andl	$0x8000, %eax
199	orl	%eax, 8(%edx)
200	fstp	%st(0)
201	ftst
202	fnstsw
203	shll	$7, %eax
204	andl	$0x8000, %eax
205	orl	%eax, 20(%edx)
206	fstp	%st(0)
207	ret	$4
208
209	/* The real part is +-Inf and the imaginary part is +-0.  So return
210	   +-Inf+-0i.  */
211	.align ALIGNARG(4)
2124:	movl	4(%esp), %eax		/* Pointer to memory for result.  */
213	fstpt	12(%eax)
214	shrl	$5, %edx
215	fstp	%st(0)
216	andl	$16, %edx
217	fldl	MOX(huge_nan_null_null,%edx,1)
218	fstpt	(%eax)
219	ret	$4
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	fldl	MO(infinity)	/* Raise invalid exception.  */
230	fmull	MO(zero)
231	fstp	%st(0)
23230:	movl	%edx, %eax
233	shrl	$5, %edx
234	shll	$4, %eax
235	andl	$16, %edx
236	andl	$32, %eax
237	orl	%eax, %edx
238	movl	4(%esp), %eax		/* Pointer to memory for result.  */
239
240	fldl	MOX(huge_nan_null_null,%edx,1)
241	fldl	MOX(huge_nan_null_null+8,%edx,1)
242	fxch
243	fstpt	(%eax)
244	fstpt	12(%eax)
245	ret	$4
246
247	/* The real part is NaN.  */
248	.align ALIGNARG(4)
24920:	fldl	MO(infinity)		/* Raise invalid exception.  */
250	fmull	MO(zero)
251	fstp	%st(0)
2522:	fstp	%st(0)
253	fstp	%st(0)
254	movl	4(%esp), %eax		/* Pointer to memory for result.  */
255	fldl	MO(huge_nan_null_null+8)
256	fld	%st(0)
257	fstpt	(%eax)
258	fstpt	12(%eax)
259	ret	$4
260
261END(__cexpl)
262weak_alias (__cexpl, cexpl)
263