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_pow) 62 fldl 12(%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 fldl 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: fldl 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 fldl 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: fldl 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 the number of bits small 226 // enough so that all are coming from the mantissa? 227 popl %eax 228 popl %edx 229 andb $1, %al 230 jz 18f // jump if not odd 231 movl %edx, %eax 232 orl %edx, %edx 233 jns 155f 234 negl %eax 235155: cmpl $0x00200000, %eax 236 ja 18f // does not fit in mantissa bits 237 // It's an odd integer. 238 shrl $31, %edx 239 fldl MOX(minf_mzero, %edx, 8) 240 ret 241 242 .align ALIGNARG(4) 24316: fcompl MO(zero) 244 addl $8, %esp 245 fnstsw 246 shrl $5, %eax 247 andl $8, %eax 248 fldl MOX(inf_zero, %eax, 1) 249 ret 250 251 .align ALIGNARG(4) 25217: shll $30, %edx // sign bit for y in right position 253 addl $8, %esp 25418: shrl $31, %edx 255 fldl MOX(inf_zero, %edx, 8) 256 ret 257 258 .align ALIGNARG(4) 259 // x is �0 26020: fstp %st(0) // y 261 testb $2, %dl 262 jz 21f // y > 0 263 264 // x is �0 and y is < 0. We must find out whether y is an odd integer. 265 testb $2, %dh 266 jz 25f 267 268 fld %st // y : y 269 fistpll (%esp) // y 270 fildll (%esp) // int(y) : y 271 fucompp // <empty> 272 fnstsw 273 sahf 274 jne 26f 275 276 // OK, the value is an integer, but is the number of bits small 277 // enough so that all are coming from the mantissa? 278 popl %eax 279 popl %edx 280 andb $1, %al 281 jz 27f // jump if not odd 282 cmpl $0xffe00000, %edx 283 jbe 27f // does not fit in mantissa bits 284 // It's an odd integer. 285 // Raise divide-by-zero exception and get minus infinity value. 286 fldl MO(one) 287 fdivl MO(zero) 288 fchs 289 ret 290 29125: fstp %st(0) 29226: addl $8, %esp 29327: // Raise divide-by-zero exception and get infinity value. 294 fldl MO(one) 295 fdivl MO(zero) 296 ret 297 298 .align ALIGNARG(4) 299 // x is �0 and y is > 0. We must find out whether y is an odd integer. 30021: testb $2, %dh 301 jz 22f 302 303 fld %st // y : y 304 fistpll (%esp) // y 305 fildll (%esp) // int(y) : y 306 fucompp // <empty> 307 fnstsw 308 sahf 309 jne 23f 310 311 // OK, the value is an integer, but is the number of bits small 312 // enough so that all are coming from the mantissa? 313 popl %eax 314 popl %edx 315 andb $1, %al 316 jz 24f // jump if not odd 317 cmpl $0xffe00000, %edx 318 jae 24f // does not fit in mantissa bits 319 // It's an odd integer. 320 fldl MO(mzero) 321 ret 322 32322: fstp %st(0) 32423: addl $8, %esp // Don't use 2 x pop 32524: fldl MO(zero) 326 ret 327 328END(__ieee754_pow) 329