1/* ix87 specific implementation of pow function. 2 Copyright (C) 1996, 1997, 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_powf) 62 flds 8(%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 flds 4(%esp) // x : y 84 85 subl $4, %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 fistpl (%esp) // y : x 103 fildl (%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 %edx 111 orl $0, %edx 112 fstp %st(0) // x 113 jns 4f // y >= 0, jump 114 fdivrl MO(one) // 1/x (now referred to as x) 115 negl %edx 1164: fldl MO(one) // 1 : x 117 fxch 118 1196: shrl $1, %edx 120 jnc 5f 121 fxch 122 fmul %st(1) // x : ST*x 123 fxch 1245: fmul %st(0), %st // x*x : ST*x 125 testl %edx, %edx 126 jnz 6b 127 fstp %st(0) // ST*x 128 ret 129 130 /* y is �NAN */ 13130: flds 4(%esp) // x : y 132 fldl MO(one) // 1.0 : x : y 133 fucomp %st(1) // x : y 134 fnstsw 135 sahf 136 je 31f 137 fxch // y : x 13831: fstp %st(1) 139 ret 140 141 .align ALIGNARG(4) 1422: /* y is a real number. */ 143 fxch // x : y 144 fldl MO(one) // 1.0 : x : y 145 fld %st(1) // x : 1.0 : x : y 146 fsub %st(1) // x-1 : 1.0 : x : y 147 fabs // |x-1| : 1.0 : x : y 148 fcompl MO(limit) // 1.0 : x : y 149 fnstsw 150 fxch // x : 1.0 : y 151 sahf 152 ja 7f 153 fsub %st(1) // x-1 : 1.0 : y 154 fyl2xp1 // log2(x) : y 155 jmp 8f 156 1577: fyl2x // log2(x) : y 1588: fmul %st(1) // y*log2(x) : y 159 fst %st(1) // y*log2(x) : y*log2(x) 160 frndint // int(y*log2(x)) : y*log2(x) 161 fsubr %st, %st(1) // int(y*log2(x)) : fract(y*log2(x)) 162 fxch // fract(y*log2(x)) : int(y*log2(x)) 163 f2xm1 // 2^fract(y*log2(x))-1 : int(y*log2(x)) 164 faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x)) 165 fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x)) 166 addl $4, %esp 167 fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x)) 168 ret 169 170 171 // pow(x,�0) = 1 172 .align ALIGNARG(4) 17311: fstp %st(0) // pop y 174 fldl MO(one) 175 ret 176 177 // y == �inf 178 .align ALIGNARG(4) 17912: fstp %st(0) // pop y 180 flds 4(%esp) // x 181 fabs 182 fcompl MO(one) // < 1, == 1, or > 1 183 fnstsw 184 andb $0x45, %ah 185 cmpb $0x45, %ah 186 je 13f // jump if x is NaN 187 188 cmpb $0x40, %ah 189 je 14f // jump if |x| == 1 190 191 shlb $1, %ah 192 xorb %ah, %dl 193 andl $2, %edx 194 fldl MOX(inf_zero, %edx, 4) 195 ret 196 197 .align ALIGNARG(4) 19814: fldl MO(one) 199 ret 200 201 .align ALIGNARG(4) 20213: flds 4(%esp) // load x == NaN 203 ret 204 205 .align ALIGNARG(4) 206 // x is �inf 20715: fstp %st(0) // y 208 testb $2, %dh 209 jz 16f // jump if x == +inf 210 211 // We must find out whether y is an odd integer. 212 fld %st // y : y 213 fistpl (%esp) // y 214 fildl (%esp) // int(y) : y 215 fucompp // <empty> 216 fnstsw 217 sahf 218 jne 17f 219 220 // OK, the value is an integer, but is the number of bits small 221 // enough so that all are coming from the mantissa? 222 popl %edx 223 testb $1, %dl 224 jz 18f // jump if not odd 225 movl %edx, %eax 226 orl %edx, %edx 227 jns 155f 228 negl %eax 229155: cmpl $0x01000000, %eax 230 ja 18f // does not fit in mantissa bits 231 // It's an odd integer. 232 shrl $31, %edx 233 fldl MOX(minf_mzero, %edx, 8) 234 ret 235 236 .align ALIGNARG(4) 23716: fcompl MO(zero) 238 addl $4, %esp 239 fnstsw 240 shrl $5, %eax 241 andl $8, %eax 242 fldl MOX(inf_zero, %eax, 1) 243 ret 244 245 .align ALIGNARG(4) 24617: shll $30, %edx // sign bit for y in right position 247 addl $4, %esp 24818: shrl $31, %edx 249 fldl MOX(inf_zero, %edx, 8) 250 ret 251 252 .align ALIGNARG(4) 253 // x is �0 25420: fstp %st(0) // y 255 testb $2, %dl 256 jz 21f // y > 0 257 258 // x is �0 and y is < 0. We must find out whether y is an odd integer. 259 testb $2, %dh 260 jz 25f 261 262 fld %st // y : y 263 fistpl (%esp) // y 264 fildl (%esp) // int(y) : y 265 fucompp // <empty> 266 fnstsw 267 sahf 268 jne 26f 269 270 // OK, the value is an integer, but is the number of bits small 271 // enough so that all are coming from the mantissa? 272 popl %edx 273 testb $1, %dl 274 jz 27f // jump if not odd 275 cmpl $0xff000000, %edx 276 jbe 27f // does not fit in mantissa bits 277 // It's an odd integer. 278 // Raise divide-by-zero exception and get minus infinity value. 279 fldl MO(one) 280 fdivl MO(zero) 281 fchs 282 ret 283 28425: fstp %st(0) 28526: addl $4, %esp 28627: // Raise divide-by-zero exception and get infinity value. 287 fldl MO(one) 288 fdivl MO(zero) 289 ret 290 291 .align ALIGNARG(4) 292 // x is �0 and y is > 0. We must find out whether y is an odd integer. 29321: testb $2, %dh 294 jz 22f 295 296 fld %st // y : y 297 fistpl (%esp) // y 298 fildl (%esp) // int(y) : y 299 fucompp // <empty> 300 fnstsw 301 sahf 302 jne 23f 303 304 // OK, the value is an integer, but is the number of bits small 305 // enough so that all are coming from the mantissa? 306 popl %edx 307 testb $1, %dl 308 jz 24f // jump if not odd 309 cmpl $0xff000000, %edx 310 jae 24f // does not fit in mantissa bits 311 // It's an odd integer. 312 fldl MO(mzero) 313 ret 314 31522: fstp %st(0) 31623: addl $4, %esp // Don't use pop. 31724: fldl MO(zero) 318 ret 319 320END(__ieee754_powf) 321