1 /* ix87 specific implementation of complex exponential function for double.
2 Copyright (C) 1997, 2005, 2012 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4 Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
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.
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.
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
26 ASM_TYPE_DIRECTIVE(huge_nan_null_null,@object)
28 .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
29 .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
33 .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
34 .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
36 .byte 0, 0, 0, 0, 0, 0, 0, 0x80
37 ASM_SIZE_DIRECTIVE(huge_nan_null_null)
39 ASM_TYPE_DIRECTIVE(twopi,@object)
41 .byte 0x35, 0xc2, 0x68, 0x21, 0xa2, 0xda, 0xf, 0xc9, 0x1, 0x40
42 .byte 0, 0, 0, 0, 0, 0
43 ASM_SIZE_DIRECTIVE(twopi)
45 ASM_TYPE_DIRECTIVE(l2e,@object)
47 .byte 0xbc, 0xf0, 0x17, 0x5c, 0x29, 0x3b, 0xaa, 0xb8, 0xff, 0x3f
48 .byte 0, 0, 0, 0, 0, 0
49 ASM_SIZE_DIRECTIVE(l2e)
51 ASM_TYPE_DIRECTIVE(one,@object)
53 ASM_SIZE_DIRECTIVE(one)
57 #define MO(op) op##@GOTOFF(%ecx)
58 #define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
61 #define MOX(op,x,f) op(,x,f)
69 fldl 16(%esp) /* y : x */
76 je 1f /* Jump if real part is +-Inf */
78 je 2f /* Jump if real part is NaN */
82 /* If the imaginary part is not finite we return NaN+i NaN, as
83 for the case when the real part is NaN. A test for +-Inf and
84 NaN would be necessary. But since we know the stack register
85 we applied `fxam' to is not empty we can simply use one test.
86 Check your FPU manual for more information. */
91 /* We have finite numbers in the real and imaginary part. Do
94 fldt MO(l2e) /* log2(e) : x : y */
95 fmulp /* x * log2(e) : y */
96 fld %st /* x * log2(e) : x * log2(e) : y */
97 frndint /* int(x * log2(e)) : x * log2(e) : y */
98 fsubr %st, %st(1) /* int(x * log2(e)) : frac(x * log2(e)) : y */
99 fxch /* frac(x * log2(e)) : int(x * log2(e)) : y */
100 f2xm1 /* 2^frac(x * log2(e))-1 : int(x * log2(e)) : y */
101 faddl MO(one) /* 2^frac(x * log2(e)) : int(x * log2(e)) : y */
102 fscale /* e^x : int(x * log2(e)) : y */
103 fst %st(1) /* e^x : e^x : y */
104 fxch %st(2) /* y : e^x : e^x */
105 fsincos /* cos(y) : sin(y) : e^x : e^x */
109 fmulp %st, %st(3) /* sin(y) : e^x : e^x * cos(y) */
110 fmulp %st, %st(1) /* e^x * sin(y) : e^x * cos(y) */
111 movl 4(%esp), %eax /* Pointer to memory for result. */
116 /* We have to reduce the argument to fsincos. */
118 7: fldt MO(twopi) /* 2*pi : y : e^x : e^x */
119 fxch /* y : 2*pi : e^x : e^x */
120 8: fprem1 /* y%(2*pi) : 2*pi : e^x : e^x */
124 fstp %st(1) /* y%(2*pi) : e^x : e^x */
125 fsincos /* cos(y) : sin(y) : e^x : e^x */
128 movl 4(%esp), %eax /* Pointer to memory for result. */
133 /* The real part is +-inf. We must make further differences. */
138 testb $0x01, %ah /* See above why 0x01 is usable here. */
142 /* The real part is +-Inf and the imaginary part is finite. */
144 cmpb $0x40, %dl /* Imaginary part == 0? */
149 fstp %st(0) /* y */ /* Drop the real part. */
150 andl $16, %edx /* This puts the sign bit of the real part
151 in bit 4. So we can use it to index a
152 small array to select 0 or Inf. */
153 fsincos /* cos(y) : sin(y) */
157 fldl MOX(huge_nan_null_null,%edx,1)
158 movl 4(%esp), %edx /* Pointer to memory for result. */
164 andl $0x80000000, %eax
170 andl $0x80000000, %eax
174 /* We must reduce the argument to fsincos. */
184 fldl MOX(huge_nan_null_null,%edx,1)
185 movl 4(%esp), %edx /* Pointer to memory for result. */
191 andl $0x80000000, %eax
197 andl $0x80000000, %eax
202 /* The real part is +-Inf and the imaginary part is +-0. So return
205 4: movl 4(%esp), %eax /* Pointer to memory for result. */
210 fldl MOX(huge_nan_null_null,%edx,1)
214 /* The real part is +-Inf, the imaginary is also is not finite. */
217 fstp %st(0) /* <empty> */
222 fldl MO(infinity) /* Raise invalid exception. */
231 movl 4(%esp), %eax /* Pointer to memory for result. */
233 fldl MOX(huge_nan_null_null,%edx,1)
234 fldl MOX(huge_nan_null_null+8,%edx,1)
240 /* The real part is NaN. */
242 20: fldl MO(infinity) /* Raise invalid exception. */
247 movl 4(%esp), %eax /* Pointer to memory for result. */
248 fldl MO(huge_nan_null_null+8)
254 weak_alias (__cexp, cexp)