1/* ix87 specific implementation of pow function. 2 Copyright (C) 1996, 1997, 1998, 1999, 2001 Free Software Foundation, Inc. 3 This file is part of the GNU C Library. 4 Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996. 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 <machine/asm.h> 22 23#ifdef __ELF__ 24 .section .rodata 25#else 26 .text 27#endif 28 29 .align ALIGNARG(4) 30 ASM_TYPE_DIRECTIVE(infinity,@object) 31inf_zero: 32infinity: 33 .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f 34 ASM_SIZE_DIRECTIVE(infinity) 35 ASM_TYPE_DIRECTIVE(zero,@object) 36zero: .double 0.0 37 ASM_SIZE_DIRECTIVE(zero) 38 ASM_TYPE_DIRECTIVE(minf_mzero,@object) 39minf_mzero: 40minfinity: 41 .byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff 42mzero: 43 .byte 0, 0, 0, 0, 0, 0, 0, 0x80 44 ASM_SIZE_DIRECTIVE(minf_mzero) 45 ASM_TYPE_DIRECTIVE(one,@object) 46one: .double 1.0 47 ASM_SIZE_DIRECTIVE(one) 48 ASM_TYPE_DIRECTIVE(limit,@object) 49limit: .double 0.29 50 ASM_SIZE_DIRECTIVE(limit) 51 52#ifdef PIC 53#define MO(op) op##@GOTOFF(%ecx) 54#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f) 55#else 56#define MO(op) op 57#define MOX(op,x,f) op(,x,f) 58#endif 59 60 .text 61ENTRY(__ieee754_powl) 62 fldt 16(%esp) // y 63 fxam 64 65#ifdef PIC 66 call 1f 671: popl %ecx 68 addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx 69#endif 70 71 fnstsw 72 movb %ah, %dl 73 andb $0x45, %ah 74 cmpb $0x40, %ah // is y == 0 ? 75 je 11f 76 77 cmpb $0x05, %ah // is y == �inf ? 78 je 12f 79 80 cmpb $0x01, %ah // is y == NaN ? 81 je 30f 82 83 fldt 4(%esp) // x : y 84 85 subl $8,%esp 86 87 fxam 88 fnstsw 89 movb %ah, %dh 90 andb $0x45, %ah 91 cmpb $0x40, %ah 92 je 20f // x is �0 93 94 cmpb $0x05, %ah 95 je 15f // x is �inf 96 97 fxch // y : x 98 99 /* First see whether `y' is a natural number. In this case we 100 can use a more precise algorithm. */ 101 fld %st // y : y : x 102 fistpll (%esp) // y : x 103 fildll (%esp) // int(y) : y : x 104 fucomp %st(1) // y : x 105 fnstsw 106 sahf 107 jne 2f 108 109 /* OK, we have an integer value for y. */ 110 popl %eax 111 popl %edx 112 orl $0, %edx 113 fstp %st(0) // x 114 jns 4f // y >= 0, jump 115 fdivrl MO(one) // 1/x (now referred to as x) 116 negl %eax 117 adcl $0, %edx 118 negl %edx 1194: fldl MO(one) // 1 : x 120 fxch 121 1226: shrdl $1, %edx, %eax 123 jnc 5f 124 fxch 125 fmul %st(1) // x : ST*x 126 fxch 1275: fmul %st(0), %st // x*x : ST*x 128 shrl $1, %edx 129 movl %eax, %ecx 130 orl %edx, %ecx 131 jnz 6b 132 fstp %st(0) // ST*x 133 ret 134 135 /* y is �NAN */ 13630: fldt 4(%esp) // x : y 137 fldl MO(one) // 1.0 : x : y 138 fucomp %st(1) // x : y 139 fnstsw 140 sahf 141 je 31f 142 fxch // y : x 14331: fstp %st(1) 144 ret 145 146 .align ALIGNARG(4) 1472: /* y is a real number. */ 148 fxch // x : y 149 fldl MO(one) // 1.0 : x : y 150 fld %st(1) // x : 1.0 : x : y 151 fsub %st(1) // x-1 : 1.0 : x : y 152 fabs // |x-1| : 1.0 : x : y 153 fcompl MO(limit) // 1.0 : x : y 154 fnstsw 155 fxch // x : 1.0 : y 156 sahf 157 ja 7f 158 fsub %st(1) // x-1 : 1.0 : y 159 fyl2xp1 // log2(x) : y 160 jmp 8f 161 1627: fyl2x // log2(x) : y 1638: fmul %st(1) // y*log2(x) : y 164 fst %st(1) // y*log2(x) : y*log2(x) 165 frndint // int(y*log2(x)) : y*log2(x) 166 fsubr %st, %st(1) // int(y*log2(x)) : fract(y*log2(x)) 167 fxch // fract(y*log2(x)) : int(y*log2(x)) 168 f2xm1 // 2^fract(y*log2(x))-1 : int(y*log2(x)) 169 faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x)) 170 fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x)) 171 addl $8, %esp 172 fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x)) 173 ret 174 175 176 // pow(x,�0) = 1 177 .align ALIGNARG(4) 17811: fstp %st(0) // pop y 179 fldl MO(one) 180 ret 181 182 // y == �inf 183 .align ALIGNARG(4) 18412: fstp %st(0) // pop y 185 fldt 4(%esp) // x 186 fabs 187 fcompl MO(one) // < 1, == 1, or > 1 188 fnstsw 189 andb $0x45, %ah 190 cmpb $0x45, %ah 191 je 13f // jump if x is NaN 192 193 cmpb $0x40, %ah 194 je 14f // jump if |x| == 1 195 196 shlb $1, %ah 197 xorb %ah, %dl 198 andl $2, %edx 199 fldl MOX(inf_zero, %edx, 4) 200 ret 201 202 .align ALIGNARG(4) 20314: fldl MO(one) 204 ret 205 206 .align ALIGNARG(4) 20713: fldt 4(%esp) // load x == NaN 208 ret 209 210 .align ALIGNARG(4) 211 // x is �inf 21215: fstp %st(0) // y 213 testb $2, %dh 214 jz 16f // jump if x == +inf 215 216 // We must find out whether y is an odd integer. 217 fld %st // y : y 218 fistpll (%esp) // y 219 fildll (%esp) // int(y) : y 220 fucompp // <empty> 221 fnstsw 222 sahf 223 jne 17f 224 225 // OK, the value is an integer, but is it odd? 226 popl %eax 227 popl %edx 228 andb $1, %al 229 jz 18f // jump if not odd 230 // It's an odd integer. 231 shrl $31, %edx 232 fldl MOX(minf_mzero, %edx, 8) 233 ret 234 235 .align ALIGNARG(4) 23616: fcompl MO(zero) 237 addl $8, %esp 238 fnstsw 239 shrl $5, %eax 240 andl $8, %eax 241 fldl MOX(inf_zero, %eax, 1) 242 ret 243 244 .align ALIGNARG(4) 24517: shll $30, %edx // sign bit for y in right position 246 addl $8, %esp 24718: shrl $31, %edx 248 fldl MOX(inf_zero, %edx, 8) 249 ret 250 251 .align ALIGNARG(4) 252 // x is �0 25320: fstp %st(0) // y 254 testb $2, %dl 255 jz 21f // y > 0 256 257 // x is �0 and y is < 0. We must find out whether y is an odd integer. 258 testb $2, %dh 259 jz 25f 260 261 fld %st // y : y 262 fistpll (%esp) // y 263 fildll (%esp) // int(y) : y 264 fucompp // <empty> 265 fnstsw 266 sahf 267 jne 26f 268 269 // OK, the value is an integer, but is it odd? 270 popl %eax 271 popl %edx 272 andb $1, %al 273 jz 27f // jump if not odd 274 // It's an odd integer. 275 // Raise divide-by-zero exception and get minus infinity value. 276 fldl MO(one) 277 fdivl MO(zero) 278 fchs 279 ret 280 28125: fstp %st(0) 28226: addl $8, %esp 28327: // Raise divide-by-zero exception and get infinity value. 284 fldl MO(one) 285 fdivl MO(zero) 286 ret 287 288 .align ALIGNARG(4) 289 // x is �0 and y is > 0. We must find out whether y is an odd integer. 29021: testb $2, %dh 291 jz 22f 292 293 fld %st // y : y 294 fistpll (%esp) // y 295 fildll (%esp) // int(y) : y 296 fucompp // <empty> 297 fnstsw 298 sahf 299 jne 23f 300 301 // OK, the value is an integer, but is it odd? 302 popl %eax 303 popl %edx 304 andb $1, %al 305 jz 24f // jump if not odd 306 // It's an odd integer. 307 fldl MO(mzero) 308 ret 309 31022: fstp %st(0) 31123: addl $8, %esp // Don't use 2 x pop 31224: fldl MO(zero) 313 ret 314 315END(__ieee754_powl) 316