From 5700f5fcc32914049a6d8c48021c75ca0f8e4fd8 Mon Sep 17 00:00:00 2001 From: Jessica Hamilton Date: Mon, 2 Jun 2014 12:57:35 +1200 Subject: [PATCH] libroot.so: update glibc's e_pow.S on x86. Fixes #9962 --- .../libroot/posix/glibc/arch/x86/e_pow.S | 213 ++++++++++++------ 1 file changed, 150 insertions(+), 63 deletions(-) diff --git a/src/system/libroot/posix/glibc/arch/x86/e_pow.S b/src/system/libroot/posix/glibc/arch/x86/e_pow.S index 997dd30e9d..76be3d138c 100644 --- a/src/system/libroot/posix/glibc/arch/x86/e_pow.S +++ b/src/system/libroot/posix/glibc/arch/x86/e_pow.S @@ -1,5 +1,13 @@ /* ix87 specific implementation of pow function. - Copyright (C) 1996, 1997, 1998, 1999, 2001 Free Software Foundation, Inc. + Copyright (C) 1996-2014 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper , 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 +/* ix87 specific implementation of pow function. + Copyright (C) 1996-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper , 1996. @@ -14,47 +22,52 @@ 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, write to the Free - Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA - 02111-1307 USA. */ + License along with the GNU C Library; if not, see + . */ #include -#ifdef __ELF__ - .section .rodata -#else - .text -#endif + .section .rodata.cst8,"aM",@progbits,8 - .align ALIGNARG(4) - ASM_TYPE_DIRECTIVE(infinity,@object) + .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) - ASM_TYPE_DIRECTIVE(zero,@object) + .type zero,@object zero: .double 0.0 ASM_SIZE_DIRECTIVE(zero) - ASM_TYPE_DIRECTIVE(minf_mzero,@object) + .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) - 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) -#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f) +# 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) +# define MO(op) op +# define MOX(op,x,f) op(,x,f) #endif .text @@ -63,9 +76,7 @@ ENTRY(__ieee754_pow) fxam #ifdef PIC - call 1f -1: popl %ecx - addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx + LOAD_PIC_REG (cx) #endif fnstsw @@ -74,7 +85,7 @@ ENTRY(__ieee754_pow) cmpb $0x40, %ah // is y == 0 ? je 11f - cmpb $0x05, %ah // is y == ±inf ? + cmpb $0x05, %ah // is y == ±inf ? je 12f cmpb $0x01, %ah // is y == NaN ? @@ -89,13 +100,24 @@ ENTRY(__ieee754_pow) movb %ah, %dh andb $0x45, %ah cmpb $0x40, %ah - je 20f // x is ±0 + je 20f // x is ±0 cmpb $0x05, %ah - je 15f // x is ±inf + 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 @@ -104,9 +126,17 @@ ENTRY(__ieee754_pow) fucomp %st(1) // y : x fnstsw sahf - jne 2f + jne 3f - /* OK, we have an integer value for y. */ + /* 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 popl %edx orl $0, %edx @@ -132,7 +162,7 @@ ENTRY(__ieee754_pow) fstp %st(0) // ST*x ret - /* y is ±NAN */ + /* y is ±NAN */ 30: fldl 4(%esp) // x : y fldl MO(one) // 1.0 : x : y fucomp %st(1) // x : y @@ -143,14 +173,26 @@ ENTRY(__ieee754_pow) 31: fstp %st(1) ret +32: addl $8, %esp + fstp %st(1) + ret + .align ALIGNARG(4) -2: /* y is a real number. */ +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 - 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 + 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 @@ -168,23 +210,53 @@ ENTRY(__ieee754_pow) 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)) + testb $2, %dh + jz 292f + // x is negative. If y is an odd integer, negate the result. + fldl 20(%esp) // y : abs(result) + fld %st // y : y : abs(result) + fabs // |y| : y : abs(result) + fcompl MO(p63) // y : abs(result) + fnstsw + sahf + jnc 291f + + // We must find out whether y is an odd integer. + fld %st // y : y : abs(result) + fistpll (%esp) // y : abs(result) + fildll (%esp) // int(y) : y : abs(result) + fucompp // abs(result) + fnstsw + sahf + jne 292f + + // OK, the value is an integer, but is it odd? + popl %eax + popl %edx + andb $1, %al + jz 290f // jump if not odd + // It's an odd integer. + fchs +290: ret +291: fstp %st(0) // abs(result) +292: addl $8, %esp ret - // pow(x,±0) = 1 + // pow(x,±0) = 1 .align ALIGNARG(4) 11: fstp %st(0) // pop y fldl MO(one) ret - // y == ±inf + // y == ±inf .align ALIGNARG(4) 12: fstp %st(0) // pop y - fldl 4(%esp) // x - fabs - fcompl MO(one) // < 1, == 1, or > 1 + 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 @@ -208,11 +280,21 @@ ENTRY(__ieee754_pow) ret .align ALIGNARG(4) - // x is ±inf + // 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 @@ -222,18 +304,11 @@ ENTRY(__ieee754_pow) sahf jne 17f - // OK, the value is an integer, but is the number of bits small - // enough so that all are coming from the mantissa? + // OK, the value is an integer. popl %eax popl %edx andb $1, %al jz 18f // jump if not odd - movl %edx, %eax - orl %edx, %edx - jns 155f - negl %eax -155: cmpl $0x00200000, %eax - ja 18f // does not fit in mantissa bits // It's an odd integer. shrl $31, %edx fldl MOX(minf_mzero, %edx, 8) @@ -256,15 +331,25 @@ ENTRY(__ieee754_pow) ret .align ALIGNARG(4) - // x is ±0 + // 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. + // 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 @@ -273,14 +358,11 @@ ENTRY(__ieee754_pow) sahf jne 26f - // OK, the value is an integer, but is the number of bits small - // enough so that all are coming from the mantissa? + // OK, the value is an integer. popl %eax popl %edx andb $1, %al jz 27f // jump if not odd - cmpl $0xffe00000, %edx - jbe 27f // does not fit in mantissa bits // It's an odd integer. // Raise divide-by-zero exception and get minus infinity value. fldl MO(one) @@ -296,10 +378,18 @@ ENTRY(__ieee754_pow) ret .align ALIGNARG(4) - // x is ±0 and y is > 0. We must find out whether y is an odd integer. + // 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 @@ -308,14 +398,11 @@ ENTRY(__ieee754_pow) sahf jne 23f - // OK, the value is an integer, but is the number of bits small - // enough so that all are coming from the mantissa? + // OK, the value is an integer. popl %eax popl %edx andb $1, %al jz 24f // jump if not odd - cmpl $0xffe00000, %edx - jae 24f // does not fit in mantissa bits // It's an odd integer. fldl MO(mzero) ret