1/* ix87 specific implementation of pow function. 2 Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004, 2007 3 Free Software Foundation, Inc. 4 This file is part of the GNU C Library. 5 Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996. 6 7 The GNU C Library is free software; you can redistribute it and/or 8 modify it under the terms of the GNU Lesser General Public 9 License as published by the Free Software Foundation; either 10 version 2.1 of the License, or (at your option) any later version. 11 12 The GNU C Library is distributed in the hope that it will be useful, 13 but WITHOUT ANY WARRANTY; without even the implied warranty of 14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 15 Lesser General Public License for more details. 16 17 You should have received a copy of the GNU Lesser General Public 18 License along with the GNU C Library; if not, write to the Free 19 Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 20 02111-1307 USA. */ 21 22#include <machine/asm.h> 23 24#ifdef __ELF__ 25 .section .rodata 26#else 27 .text 28#endif 29 30 .align ALIGNARG(4) 31 ASM_TYPE_DIRECTIVE(infinity,@object) 32inf_zero: 33infinity: 34 .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f 35 ASM_SIZE_DIRECTIVE(infinity) 36 ASM_TYPE_DIRECTIVE(zero,@object) 37zero: .double 0.0 38 ASM_SIZE_DIRECTIVE(zero) 39 ASM_TYPE_DIRECTIVE(minf_mzero,@object) 40minf_mzero: 41minfinity: 42 .byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff 43mzero: 44 .byte 0, 0, 0, 0, 0, 0, 0, 0x80 45 ASM_SIZE_DIRECTIVE(minf_mzero) 46 ASM_TYPE_DIRECTIVE(one,@object) 47one: .double 1.0 48 ASM_SIZE_DIRECTIVE(one) 49 ASM_TYPE_DIRECTIVE(limit,@object) 50limit: .double 0.29 51 ASM_SIZE_DIRECTIVE(limit) 52 ASM_TYPE_DIRECTIVE(p63,@object) 53p63: 54 .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43 55 ASM_SIZE_DIRECTIVE(p63) 56 57#ifdef PIC 58#define MO(op) op##(%rip) 59#else 60#define MO(op) op 61#endif 62 63 .text 64ENTRY(__ieee754_powl) 65 fldt 24(%rsp) // y 66 fxam 67 68 69 fnstsw 70 movb %ah, %dl 71 andb $0x45, %ah 72 cmpb $0x40, %ah // is y == 0 ? 73 je 11f 74 75 cmpb $0x05, %ah // is y == �inf ? 76 je 12f 77 78 cmpb $0x01, %ah // is y == NaN ? 79 je 30f 80 81 fldt 8(%rsp) // x : y 82 83 fxam 84 fnstsw 85 movb %ah, %dh 86 andb $0x45, %ah 87 cmpb $0x40, %ah 88 je 20f // x is �0 89 90 cmpb $0x05, %ah 91 je 15f // x is �inf 92 93 fxch // y : x 94 95 /* fistpll raises invalid exception for |y| >= 1L<<63. */ 96 fldl MO(p63) // 1L<<63 : y : x 97 fld %st(1) // y : 1L<<63 : y : x 98 fabs // |y| : 1L<<63 : y : x 99 fcomip %st(1), %st // 1L<<63 : y : x 100 fstp %st(0) // y : x 101 jnc 2f 102 103 /* First see whether `y' is a natural number. In this case we 104 can use a more precise algorithm. */ 105 fld %st // y : y : x 106 fistpll -8(%rsp) // y : x 107 fildll -8(%rsp) // int(y) : y : x 108 fucomip %st(1),%st // y : x 109 jne 2f 110 111 /* OK, we have an integer value for y. */ 112 mov -8(%rsp),%eax 113 mov -4(%rsp),%edx 114 orl $0, %edx 115 fstp %st(0) // x 116 jns 4f // y >= 0, jump 117 fdivrl MO(one) // 1/x (now referred to as x) 118 negl %eax 119 adcl $0, %edx 120 negl %edx 1214: fldl MO(one) // 1 : x 122 fxch 123 1246: shrdl $1, %edx, %eax 125 jnc 5f 126 fxch 127 fmul %st(1) // x : ST*x 128 fxch 1295: fmul %st(0), %st // x*x : ST*x 130 shrl $1, %edx 131 movl %eax, %ecx 132 orl %edx, %ecx 133 jnz 6b 134 fstp %st(0) // ST*x 135 ret 136 137 /* y is �NAN */ 13830: fldt 8(%rsp) // x : y 139 fldl MO(one) // 1.0 : x : y 140 fucomip %st(1),%st // x : y 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 fldl MO(limit) // 0.29 : 1.0 : x : y 151 fld %st(2) // x : 0.29 : 1.0 : x : y 152 fsub %st(2) // x-1 : 0.29 : 1.0 : x : y 153 fabs // |x-1| : 0.29 : 1.0 : x : y 154 fucompp // 1.0 : x : y 155 fnstsw 156 fxch // x : 1.0 : y 157 test $4500,%eax 158 jz 7f 159 fsub %st(1) // x-1 : 1.0 : y 160 fyl2xp1 // log2(x) : y 161 jmp 8f 162 1637: fyl2x // log2(x) : y 1648: fmul %st(1) // y*log2(x) : y 165 fxam 166 fnstsw 167 andb $0x45, %ah 168 cmpb $0x05, %ah // is y*log2(x) == �inf ? 169 je 28f 170 fst %st(1) // y*log2(x) : y*log2(x) 171 frndint // int(y*log2(x)) : y*log2(x) 172 fsubr %st, %st(1) // int(y*log2(x)) : fract(y*log2(x)) 173 fxch // fract(y*log2(x)) : int(y*log2(x)) 174 f2xm1 // 2^fract(y*log2(x))-1 : int(y*log2(x)) 175 faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x)) 176 fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x)) 177 fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x)) 178 ret 179 18028: fstp %st(1) // y*log2(x) 181 fldl MO(one) // 1 : y*log2(x) 182 fscale // 2^(y*log2(x)) : y*log2(x) 183 fstp %st(1) // 2^(y*log2(x)) 184 ret 185 186 // pow(x,�0) = 1 187 .align ALIGNARG(4) 18811: fstp %st(0) // pop y 189 fldl MO(one) 190 ret 191 192 // y == �inf 193 .align ALIGNARG(4) 19412: fstp %st(0) // pop y 195 fldl MO(one) // 1 196 fldt 8(%rsp) // x : 1 197 fabs // abs(x) : 1 198 fucompp // < 1, == 1, or > 1 199 fnstsw 200 andb $0x45, %ah 201 cmpb $0x45, %ah 202 je 13f // jump if x is NaN 203 204 cmpb $0x40, %ah 205 je 14f // jump if |x| == 1 206 207 shlb $1, %ah 208 xorb %ah, %dl 209 andl $2, %edx 210#ifdef PIC 211 lea inf_zero(%rip),%rcx 212 fldl (%rcx, %rdx, 4) 213#else 214 fldl inf_zero(,%rdx, 4) 215#endif 216 ret 217 218 .align ALIGNARG(4) 21914: fldl MO(one) 220 ret 221 222 .align ALIGNARG(4) 22313: fldt 8(%rsp) // load x == NaN 224 ret 225 226 .align ALIGNARG(4) 227 // x is �inf 22815: fstp %st(0) // y 229 testb $2, %dh 230 jz 16f // jump if x == +inf 231 232 // We must find out whether y is an odd integer. 233 fld %st // y : y 234 fistpll -8(%rsp) // y 235 fildll -8(%rsp) // int(y) : y 236 fucomip %st(1),%st 237 ffreep %st // <empty> 238 jne 17f 239 240 // OK, the value is an integer, but is it odd? 241 mov -8(%rsp), %eax 242 mov -4(%rsp), %edx 243 andb $1, %al 244 jz 18f // jump if not odd 245 // It's an odd integer. 246 shrl $31, %edx 247#ifdef PIC 248 lea minf_mzero(%rip),%rcx 249 fldl (%rcx, %rdx, 8) 250#else 251 fldl minf_mzero(,%rdx, 8) 252#endif 253 ret 254 255 .align ALIGNARG(4) 25616: fcompl MO(zero) 257 fnstsw 258 shrl $5, %eax 259 andl $8, %eax 260#ifdef PIC 261 lea inf_zero(%rip),%rcx 262 fldl (%rcx, %rax, 1) 263#else 264 fldl inf_zero(,%rax, 1) 265#endif 266 ret 267 268 .align ALIGNARG(4) 26917: shll $30, %edx // sign bit for y in right position 27018: shrl $31, %edx 271#ifdef PIC 272 lea inf_zero(%rip),%rcx 273 fldl (%rcx, %rdx, 8) 274#else 275 fldl inf_zero(,%rdx, 8) 276#endif 277 ret 278 279 .align ALIGNARG(4) 280 // x is �0 28120: fstp %st(0) // y 282 testb $2, %dl 283 jz 21f // y > 0 284 285 // x is �0 and y is < 0. We must find out whether y is an odd integer. 286 testb $2, %dh 287 jz 25f 288 289 fld %st // y : y 290 fistpll -8(%rsp) // y 291 fildll -8(%rsp) // int(y) : y 292 fucomip %st(1),%st 293 ffreep %st // <empty> 294 jne 26f 295 296 // OK, the value is an integer, but is it odd? 297 mov -8(%rsp),%eax 298 mov -4(%rsp),%edx 299 andb $1, %al 300 jz 27f // jump if not odd 301 // It's an odd integer. 302 // Raise divide-by-zero exception and get minus infinity value. 303 fldl MO(one) 304 fdivl MO(zero) 305 fchs 306 ret 307 30825: fstp %st(0) 30926: 31027: // Raise divide-by-zero exception and get infinity value. 311 fldl MO(one) 312 fdivl MO(zero) 313 ret 314 315 .align ALIGNARG(4) 316 // x is �0 and y is > 0. We must find out whether y is an odd integer. 31721: testb $2, %dh 318 jz 22f 319 320 fld %st // y : y 321 fistpll -8(%rsp) // y 322 fildll -8(%rsp) // int(y) : y 323 fucomip %st(1),%st 324 ffreep %st // <empty> 325 jne 23f 326 327 // OK, the value is an integer, but is it odd? 328 mov -8(%rsp),%eax 329 mov -4(%rsp),%edx 330 andb $1, %al 331 jz 24f // jump if not odd 332 // It's an odd integer. 333 fldl MO(mzero) 334 ret 335 33622: fstp %st(0) 33723: 33824: fldl MO(zero) 339 ret 340 341END(__ieee754_powl) 342