about summary refs log tree commit diff
path: root/REORG.TODO/sysdeps/i386/fpu/e_pow.S
diff options
context:
space:
mode:
authorZack Weinberg <zackw@panix.com>2017-06-08 15:39:03 -0400
committerZack Weinberg <zackw@panix.com>2017-06-08 15:39:03 -0400
commit5046dbb4a7eba5eccfd258f92f4735c9ffc8d069 (patch)
tree4470480d904b65cf14ca524f96f79eca818c3eaf /REORG.TODO/sysdeps/i386/fpu/e_pow.S
parent199fc19d3aaaf57944ef036e15904febe877fc93 (diff)
downloadglibc-5046dbb4a7eba5eccfd258f92f4735c9ffc8d069.tar.gz
glibc-5046dbb4a7eba5eccfd258f92f4735c9ffc8d069.tar.xz
glibc-5046dbb4a7eba5eccfd258f92f4735c9ffc8d069.zip
Prepare for radical source tree reorganization. zack/build-layout-experiment
All top-level files and directories are moved into a temporary storage
directory, REORG.TODO, except for files that will certainly still
exist in their current form at top level when we're done (COPYING,
COPYING.LIB, LICENSES, NEWS, README), all old ChangeLog files (which
are moved to the new directory OldChangeLogs, instead), and the
generated file INSTALL (which is just deleted; in the new order, there
will be no generated files checked into version control).
Diffstat (limited to 'REORG.TODO/sysdeps/i386/fpu/e_pow.S')
-rw-r--r--REORG.TODO/sysdeps/i386/fpu/e_pow.S456
1 files changed, 456 insertions, 0 deletions
diff --git a/REORG.TODO/sysdeps/i386/fpu/e_pow.S b/REORG.TODO/sysdeps/i386/fpu/e_pow.S
new file mode 100644
index 0000000000..2edb9a9fbc
--- /dev/null
+++ b/REORG.TODO/sysdeps/i386/fpu/e_pow.S
@@ -0,0 +1,456 @@
+/* ix87 specific implementation of pow function.
+   Copyright (C) 1996-2017 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 Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 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
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <machine/asm.h>
+#include <i386-math-asm.h>
+
+	.section .rodata.cst8,"aM",@progbits,8
+
+	.p2align 3
+	.type one,@object
+one:	.double 1.0
+	ASM_SIZE_DIRECTIVE(one)
+	.type limit,@object
+limit:	.double 0.29
+	ASM_SIZE_DIRECTIVE(limit)
+	.type p63,@object
+p63:	.byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
+	ASM_SIZE_DIRECTIVE(p63)
+	.type p10,@object
+p10:	.byte 0, 0, 0, 0, 0, 0, 0x90, 0x40
+	ASM_SIZE_DIRECTIVE(p10)
+
+	.section .rodata.cst16,"aM",@progbits,16
+
+	.p2align 3
+	.type infinity,@object
+inf_zero:
+infinity:
+	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
+	ASM_SIZE_DIRECTIVE(infinity)
+	.type zero,@object
+zero:	.double 0.0
+	ASM_SIZE_DIRECTIVE(zero)
+	.type minf_mzero,@object
+minf_mzero:
+minfinity:
+	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff
+mzero:
+	.byte 0, 0, 0, 0, 0, 0, 0, 0x80
+	ASM_SIZE_DIRECTIVE(minf_mzero)
+DEFINE_DBL_MIN
+
+#ifdef PIC
+# define MO(op) op##@GOTOFF(%ecx)
+# define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
+#else
+# define MO(op) op
+# define MOX(op,x,f) op(,x,f)
+#endif
+
+	.text
+ENTRY(__ieee754_pow)
+	fldl	12(%esp)	// y
+	fxam
+
+#ifdef	PIC
+	LOAD_PIC_REG (cx)
+#endif
+
+	fnstsw
+	movb	%ah, %dl
+	andb	$0x45, %ah
+	cmpb	$0x40, %ah	// is y == 0 ?
+	je	11f
+
+	cmpb	$0x05, %ah	// is y == ±inf ?
+	je	12f
+
+	cmpb	$0x01, %ah	// is y == NaN ?
+	je	30f
+
+	fldl	4(%esp)		// x : y
+
+	subl	$8,%esp
+	cfi_adjust_cfa_offset (8)
+
+	fxam
+	fnstsw
+	movb	%ah, %dh
+	andb	$0x45, %ah
+	cmpb	$0x40, %ah
+	je	20f		// x is ±0
+
+	cmpb	$0x05, %ah
+	je	15f		// x is ±inf
+
+	cmpb	$0x01, %ah
+	je	32f		// x is NaN
+
+	fxch			// y : x
+
+	/* fistpll raises invalid exception for |y| >= 1L<<63.  */
+	fld	%st		// y : y : x
+	fabs			// |y| : y : x
+	fcompl	MO(p63)		// y : x
+	fnstsw
+	sahf
+	jnc	2f
+
+	/* 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	3f
+
+	/* OK, we have an integer value for y.  If large enough that
+	   errors may propagate out of the 11 bits excess precision, use
+	   the algorithm for real exponent instead.  */
+	fld	%st		// y : y : x
+	fabs			// |y| : y : x
+	fcompl	MO(p10)		// y : x
+	fnstsw
+	sahf
+	jnc	2f
+	popl	%eax
+	cfi_adjust_cfa_offset (-4)
+	popl	%edx
+	cfi_adjust_cfa_offset (-4)
+	orl	$0, %edx
+	fstp	%st(0)		// x
+	jns	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
+
+	/* If y is even, take the absolute value of x.  Otherwise,
+	   ensure all intermediate values that might overflow have the
+	   sign of x.  */
+	testb	$1, %al
+	jnz	6f
+	fabs
+
+6:	shrdl	$1, %edx, %eax
+	jnc	5f
+	fxch
+	fabs
+	fmul	%st(1)		// x : ST*x
+	fxch
+5:	fld	%st		// x : x : ST*x
+	fabs			// |x| : x : ST*x
+	fmulp			// |x|*x : ST*x
+	shrl	$1, %edx
+	movl	%eax, %ecx
+	orl	%edx, %ecx
+	jnz	6b
+	fstp	%st(0)		// ST*x
+#ifdef	PIC
+	LOAD_PIC_REG (cx)
+#endif
+	DBL_NARROW_EVAL_UFLOW_NONNAN
+	ret
+
+	/* y is ±NAN */
+30:	fldl	4(%esp)		// x : y
+	fldl	MO(one)		// 1.0 : x : y
+	fucomp	%st(1)		// x : y
+	fnstsw
+	sahf
+	je	31f
+	fxch			// y : x
+31:	fstp	%st(1)
+	ret
+
+	cfi_adjust_cfa_offset (8)
+32:	addl	$8, %esp
+	cfi_adjust_cfa_offset (-8)
+	fstp	%st(1)
+	ret
+
+	cfi_adjust_cfa_offset (8)
+	.align ALIGNARG(4)
+2:	// y is a large integer (absolute value at least 1L<<10), but
+	// may be odd unless at least 1L<<64.  So it may be necessary
+	// to adjust the sign of a negative result afterwards.
+	fxch			// x : y
+	fabs			// |x| : y
+	fxch			// y : x
+	.align ALIGNARG(4)
+3:	/* y is a real number.  */
+	fxch			// x : y
+	fldl	MO(one)		// 1.0 : x : y
+	fldl	MO(limit)	// 0.29 : 1.0 : x : y
+	fld	%st(2)		// x : 0.29 : 1.0 : x : y
+	fsub	%st(2)		// x-1 : 0.29 : 1.0 : x : y
+	fabs			// |x-1| : 0.29 : 1.0 : x : y
+	fucompp			// 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))
+
+	// Before scaling, we must negate if x is negative and y is an
+	// odd integer.
+	testb	$2, %dh
+	jz	291f
+	// x is negative.  If y is an odd integer, negate the result.
+	fldl	20(%esp)	// y : 2^fract(y*log2(x)) : int(y*log2(x))
+	fld	%st		// y : y : 2^fract(y*log2(x)) : int(y*log2(x))
+	fabs			// |y| : y : 2^fract(y*log2(x)) : int(y*log2(x))
+	fcompl	MO(p63)		// y : 2^fract(y*log2(x)) : int(y*log2(x))
+	fnstsw
+	sahf
+	jnc	290f
+
+	// We must find out whether y is an odd integer.
+	fld	%st		// y : y : 2^fract(y*log2(x)) : int(y*log2(x))
+	fistpll	(%esp)		// y : 2^fract(y*log2(x)) : int(y*log2(x))
+	fildll	(%esp)		// int(y) : y : 2^fract(y*log2(x)) : int(y*log2(x))
+	fucompp			// 2^fract(y*log2(x)) : int(y*log2(x))
+	fnstsw
+	sahf
+	jne	291f
+
+	// OK, the value is an integer, but is it odd?
+	popl	%eax
+	cfi_adjust_cfa_offset (-4)
+	popl	%edx
+	cfi_adjust_cfa_offset (-4)
+	andb	$1, %al
+	jz	292f		// jump if not odd
+	// It's an odd integer.
+	fchs
+	jmp	292f
+
+	cfi_adjust_cfa_offset (8)
+290:	fstp	%st(0)		// 2^fract(y*log2(x)) : int(y*log2(x))
+291:	addl	$8, %esp
+	cfi_adjust_cfa_offset (-8)
+292:	fscale			// +/- 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
+	fstp	%st(1)		// +/- 2^fract(y*log2(x))*2^int(y*log2(x))
+	DBL_NARROW_EVAL_UFLOW_NONNAN
+	ret
+
+
+	// pow(x,±0) = 1
+	.align ALIGNARG(4)
+11:	fstp	%st(0)		// pop y
+	fldl	MO(one)
+	ret
+
+	// y == ±inf
+	.align ALIGNARG(4)
+12:	fstp	%st(0)		// pop y
+	fldl	MO(one)		// 1
+	fldl	4(%esp)		// x : 1
+	fabs			// abs(x) : 1
+	fucompp			// < 1, == 1, or > 1
+	fnstsw
+	andb	$0x45, %ah
+	cmpb	$0x45, %ah
+	je	13f		// jump if x is NaN
+
+	cmpb	$0x40, %ah
+	je	14f		// jump if |x| == 1
+
+	shlb	$1, %ah
+	xorb	%ah, %dl
+	andl	$2, %edx
+	fldl	MOX(inf_zero, %edx, 4)
+	ret
+
+	.align ALIGNARG(4)
+14:	fldl	MO(one)
+	ret
+
+	.align ALIGNARG(4)
+13:	fldl	4(%esp)		// load x == NaN
+	ret
+
+	cfi_adjust_cfa_offset (8)
+	.align ALIGNARG(4)
+	// x is ±inf
+15:	fstp	%st(0)		// y
+	testb	$2, %dh
+	jz	16f		// jump if x == +inf
+
+	// fistpll raises invalid exception for |y| >= 1L<<63, so test
+	// that (in which case y is certainly even) before testing
+	// whether y is odd.
+	fld	%st		// y : y
+	fabs			// |y| : y
+	fcompl	MO(p63)		// y
+	fnstsw
+	sahf
+	jnc	16f
+
+	// We must find out whether y is an odd integer.
+	fld	%st		// y : y
+	fistpll	(%esp)		// y
+	fildll	(%esp)		// int(y) : y
+	fucompp			// <empty>
+	fnstsw
+	sahf
+	jne	17f
+
+	// OK, the value is an integer.
+	popl	%eax
+	cfi_adjust_cfa_offset (-4)
+	popl	%edx
+	cfi_adjust_cfa_offset (-4)
+	andb	$1, %al
+	jz	18f		// jump if not odd
+	// It's an odd integer.
+	shrl	$31, %edx
+	fldl	MOX(minf_mzero, %edx, 8)
+	ret
+
+	cfi_adjust_cfa_offset (8)
+	.align ALIGNARG(4)
+16:	fcompl	MO(zero)
+	addl	$8, %esp
+	cfi_adjust_cfa_offset (-8)
+	fnstsw
+	shrl	$5, %eax
+	andl	$8, %eax
+	fldl	MOX(inf_zero, %eax, 1)
+	ret
+
+	cfi_adjust_cfa_offset (8)
+	.align ALIGNARG(4)
+17:	shll	$30, %edx	// sign bit for y in right position
+	addl	$8, %esp
+	cfi_adjust_cfa_offset (-8)
+18:	shrl	$31, %edx
+	fldl	MOX(inf_zero, %edx, 8)
+	ret
+
+	cfi_adjust_cfa_offset (8)
+	.align ALIGNARG(4)
+	// x is ±0
+20:	fstp	%st(0)		// y
+	testb	$2, %dl
+	jz	21f		// y > 0
+
+	// x is ±0 and y is < 0.  We must find out whether y is an odd integer.
+	testb	$2, %dh
+	jz	25f
+
+	// fistpll raises invalid exception for |y| >= 1L<<63, so test
+	// that (in which case y is certainly even) before testing
+	// whether y is odd.
+	fld	%st		// y : y
+	fabs			// |y| : y
+	fcompl	MO(p63)		// y
+	fnstsw
+	sahf
+	jnc	25f
+
+	fld	%st		// y : y
+	fistpll	(%esp)		// y
+	fildll	(%esp)		// int(y) : y
+	fucompp			// <empty>
+	fnstsw
+	sahf
+	jne	26f
+
+	// OK, the value is an integer.
+	popl	%eax
+	cfi_adjust_cfa_offset (-4)
+	popl	%edx
+	cfi_adjust_cfa_offset (-4)
+	andb	$1, %al
+	jz	27f		// jump if not odd
+	// It's an odd integer.
+	// Raise divide-by-zero exception and get minus infinity value.
+	fldl	MO(one)
+	fdivl	MO(zero)
+	fchs
+	ret
+
+	cfi_adjust_cfa_offset (8)
+25:	fstp	%st(0)
+26:	addl	$8, %esp
+	cfi_adjust_cfa_offset (-8)
+27:	// Raise divide-by-zero exception and get infinity value.
+	fldl	MO(one)
+	fdivl	MO(zero)
+	ret
+
+	cfi_adjust_cfa_offset (8)
+	.align ALIGNARG(4)
+	// x is ±0 and y is > 0.  We must find out whether y is an odd integer.
+21:	testb	$2, %dh
+	jz	22f
+
+	// fistpll raises invalid exception for |y| >= 1L<<63, so test
+	// that (in which case y is certainly even) before testing
+	// whether y is odd.
+	fcoml	MO(p63)		// y
+	fnstsw
+	sahf
+	jnc	22f
+
+	fld	%st		// y : y
+	fistpll	(%esp)		// y
+	fildll	(%esp)		// int(y) : y
+	fucompp			// <empty>
+	fnstsw
+	sahf
+	jne	23f
+
+	// OK, the value is an integer.
+	popl	%eax
+	cfi_adjust_cfa_offset (-4)
+	popl	%edx
+	cfi_adjust_cfa_offset (-4)
+	andb	$1, %al
+	jz	24f		// jump if not odd
+	// It's an odd integer.
+	fldl	MO(mzero)
+	ret
+
+	cfi_adjust_cfa_offset (8)
+22:	fstp	%st(0)
+23:	addl	$8, %esp	// Don't use 2 x pop
+	cfi_adjust_cfa_offset (-8)
+24:	fldl	MO(zero)
+	ret
+
+END(__ieee754_pow)
+strong_alias (__ieee754_pow, __pow_finite)