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