]>
Commit | Line | Data |
---|---|---|
993b3242 UD |
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 Library General Public License as | |
8 | published by the Free Software Foundation; either version 2 of the | |
9 | 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 | Library General Public License for more details. | |
15 | ||
16 | You should have received a copy of the GNU Library General Public | |
17 | License along with the GNU C Library; see the file COPYING.LIB. If not, | |
18 | write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, | |
19 | Boston, MA 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) | |
30 | huge_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 | |
779ae82e UD |
34 | zero: .double 0.0 |
35 | infinity: | |
993b3242 UD |
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) | |
43 | twopi: | |
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) | |
49 | l2e: | |
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) | |
55 | one: .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 | |
68 | ENTRY(__cexp) | |
69 | fldl 8(%esp) /* x */ | |
70 | fxam | |
71 | fnstsw | |
72 | fldl 16(%esp) /* y : x */ | |
73 | #ifdef PIC | |
74 | call 1f | |
75 | 1: 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 | |
779ae82e | 94 | je 20f |
993b3242 UD |
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 | fstpl 8(%eax) | |
118 | fstpl (%eax) | |
119 | ret $4 | |
120 | ||
121 | /* We have to reduce the argument to fsincos. */ | |
122 | .align ALIGNARG(4) | |
123 | 7: fldt MO(twopi) /* 2*pi : y : e^x : e^x */ | |
124 | fxch /* y : 2*pi : e^x : e^x */ | |
125 | 8: 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 | fstpl 8(%eax) | |
135 | fstpl (%eax) | |
136 | ret $4 | |
137 | ||
138 | /* The real part is +-inf. We must make further differences. */ | |
139 | .align ALIGNARG(4) | |
140 | 1: fxam /* y : x */ | |
141 | fnstsw | |
142 | movb %ah, %dl | |
779ae82e UD |
143 | testb $0x01, %ah /* See above why 0x01 is usable here. */ |
144 | jne 3f | |
993b3242 UD |
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 | fstl 8(%edx) | |
165 | fstpl (%edx) | |
166 | ftst | |
167 | fnstsw | |
168 | shll $23, %eax | |
169 | andl $0x80000000, %eax | |
170 | orl %eax, 4(%edx) | |
171 | fstp %st(0) | |
172 | ftst | |
173 | fnstsw | |
174 | shll $23, %eax | |
175 | andl $0x80000000, %eax | |
176 | orl %eax, 12(%edx) | |
177 | fstp %st(0) | |
178 | ret $4 | |
179 | /* We must reduce the argument to fsincos. */ | |
180 | .align ALIGNARG(4) | |
181 | 5: fldt MO(twopi) | |
182 | fxch | |
183 | 6: fprem1 | |
184 | fnstsw | |
185 | testl $0x400, %eax | |
186 | jnz 6b | |
187 | fstp %st(1) | |
188 | fsincos | |
189 | fldl MOX(huge_nan_null_null,%edx,1) | |
190 | movl 4(%esp), %edx /* Pointer to memory for result. */ | |
191 | fstl 8(%edx) | |
192 | fstpl (%edx) | |
193 | ftst | |
194 | fnstsw | |
195 | shll $23, %eax | |
196 | andl $0x80000000, %eax | |
197 | orl %eax, 4(%edx) | |
198 | fstp %st(0) | |
199 | ftst | |
200 | fnstsw | |
201 | shll $23, %eax | |
202 | andl $0x80000000, %eax | |
203 | orl %eax, 12(%edx) | |
204 | fstp %st(0) | |
205 | ret $4 | |
206 | ||
207 | /* The real part is +-Inf and the imaginary part is +-0. So return | |
208 | +-Inf+-0i. */ | |
209 | .align ALIGNARG(4) | |
210 | 4: movl 4(%esp), %eax /* Pointer to memory for result. */ | |
211 | fstpl 8(%eax) | |
212 | shrl $5, %edx | |
213 | fstp %st(0) | |
214 | andl $16, %edx | |
215 | fldl MOX(huge_nan_null_null,%edx,1) | |
216 | fstpl (%eax) | |
217 | ret $4 | |
218 | ||
219 | /* The real part is +-Inf, the imaginary is also is not finite. */ | |
220 | .align ALIGNARG(4) | |
221 | 3: fstp %st(0) | |
222 | fstp %st(0) /* <empty> */ | |
779ae82e UD |
223 | andb $0x45, %ah |
224 | andb $0x47, %dh | |
225 | xorb %dh, %ah | |
226 | jnz 30f | |
227 | fldl MO(infinity) /* Raise invalid exception. */ | |
228 | fmull MO(zero) | |
229 | fstp %st(0) | |
230 | 30: movl %edx, %eax | |
993b3242 UD |
231 | shrl $5, %edx |
232 | shll $4, %eax | |
233 | andl $16, %edx | |
234 | andl $32, %eax | |
235 | orl %eax, %edx | |
236 | movl 4(%esp), %eax /* Pointer to memory for result. */ | |
237 | ||
238 | fldl MOX(huge_nan_null_null,%edx,1) | |
239 | fldl MOX(huge_nan_null_null+8,%edx,1) | |
240 | fstpl 8(%eax) | |
241 | fstpl (%eax) | |
242 | ret $4 | |
243 | ||
244 | /* The real part is NaN. */ | |
245 | .align ALIGNARG(4) | |
779ae82e UD |
246 | 20: fldl MO(infinity) /* Raise invalid exception. */ |
247 | fmull MO(zero) | |
248 | fstp %st(0) | |
993b3242 UD |
249 | 2: fstp %st(0) |
250 | fstp %st(0) | |
251 | movl 4(%esp), %eax /* Pointer to memory for result. */ | |
252 | fldl MO(huge_nan_null_null+8) | |
253 | fstl (%eax) | |
254 | fstpl 8(%eax) | |
255 | ret $4 | |
256 | ||
257 | END(__cexp) | |
258 | weak_alias (__cexp, cexp) |