aboutsummaryrefslogtreecommitdiff
path: root/sysdeps/libm-i387/e_powf.S
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/libm-i387/e_powf.S')
-rw-r--r--sysdeps/libm-i387/e_powf.S120
1 files changed, 120 insertions, 0 deletions
diff --git a/sysdeps/libm-i387/e_powf.S b/sysdeps/libm-i387/e_powf.S
new file mode 100644
index 0000000000..9d6bc510b6
--- /dev/null
+++ b/sysdeps/libm-i387/e_powf.S
@@ -0,0 +1,120 @@
+/* ix87 specific implementation of pow function.
+ Copyright (C) 1996 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Library General Public License as
+ published by the Free Software Foundation; either version 2 of the
+ License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Library General Public License for more details.
+
+ You should have received a copy of the GNU Library General Public
+ License along with the GNU C Library; see the file COPYING.LIB. If not,
+ write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+ Boston, MA 02111-1307, USA. */
+
+#include <machine/asm.h>
+
+#ifdef __ELF__
+ .section .rodata
+#else
+ .text
+#endif
+
+ .align ALIGNARG(4)
+ ASM_TYPE_DIRECTIVE(one,@object)
+one: .double 1.0
+ ASM_SIZE_DIRECTIVE(one)
+ ASM_TYPE_DIRECTIVE(limit,@object)
+limit: .double 0.29
+ ASM_SIZE_DIRECTIVE(limit)
+
+#ifdef PIC
+#define MO(op) op##@GOTOFF(%ecx)
+#else
+#define MO(op) op
+#endif
+
+ .text
+ENTRY(__ieee754_powf)
+ flds 4(%esp) // x
+ flds 8(%esp) // y : x
+
+#ifdef PIC
+ call 1f
+1: popl %ecx
+ addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx
+#endif
+ subl $8,%esp
+
+ /* First see whether `y' is a natural number. In this case we
+ can use a more precise algorithm. */
+ fld %st // y : y : x
+ fistpll (%esp) // y : x
+ fildll (%esp) // int(y) : y : x
+ fucomp %st(1) // y : x
+ fnstsw
+ sahf
+ jne 2f
+
+ /* OK, we have an integer value for y. */
+ ftst // y : x
+ fstp %st(0) // x
+ fnstsw
+ sahf
+ popl %eax
+ popl %edx
+ jnc 4f // y >= 0, jump
+ fdivrl MO(one) // 1/x (now referred to as x)
+ negl %eax
+ adcl $0, %edx
+ negl %edx
+4: fldl MO(one) // 1 : x
+ fxch
+
+6: shrdl $1, %edx, %eax
+ jnc 5f
+ fxch
+ fmul %st(1) // x : ST*x
+ fxch
+5: fmul %st(0), %st // x*x : ST*x
+ movl %eax, %ecx
+ orl %edx, %ecx
+ jnz 6b
+ fstp %st(0) // ST*x
+ ret
+
+ .align ALIGNARG(4)
+2: /* y is a real number. */
+ fxch // x : y
+ fldl MO(one) // 1.0 : x : y
+ fld %st(1) // x : 1.0 : x : y
+ fsub %st(1) // x-1 : 1.0 : x : y
+ fabs // |x-1| : 1.0 : x : y
+ fcompl MO(limit) // 1.0 : x : y
+ fnstsw
+ fxch // x : 1.0 : y
+ sahf
+ ja 7f
+ fsub %st(1) // x-1 : 1.0 : y
+ fyl2xp1 // log2(x) : y
+ jmp 8f
+
+7: fyl2x // log2(x) : y
+8: fmul %st(1) // y*log2(x) : y
+ fst %st(1) // y*log2(x) : y*log2(x)
+ frndint // int(y*log2(x)) : y*log2(x)
+ fsubr %st, %st(1) // int(y*log2(x)) : fract(y*log2(x))
+ fxch // fract(y*log2(x)) : int(y*log2(x))
+ f2xm1 // 2^fract(y*log2(x))-1 : int(y*log2(x))
+ faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x))
+ fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
+ addl $8, %esp
+ fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x))
+ ret
+END(__ieee754_powf)