diff --git a/src/system/libroot/posix/glibc/arch/x86_64/Jamfile b/src/system/libroot/posix/glibc/arch/x86_64/Jamfile index 15129e35ff..80f50cee88 100644 --- a/src/system/libroot/posix/glibc/arch/x86_64/Jamfile +++ b/src/system/libroot/posix/glibc/arch/x86_64/Jamfile @@ -23,7 +23,7 @@ SubDirCcFlags -D_GNU_SOURCE -D_IEEE_LIBM -DPIC ; SubDirAsFlags -DPIC ; local genericSources = - cmp.c dbl2mpn.c divrem.c + cmp.c dbl2mpn.c divrem.c memrchr.c mpn2dbl.c mpn2flt.c mpn2ldbl.c mul.c mul_n.c @@ -90,7 +90,7 @@ for architectureObject in [ MultiArchSubDirSetup x86_64 ] { local architecture = $(TARGET_PACKAGING_ARCH) ; MergeObject <$(architecture)>posix_gnu_arch_$(TARGET_ARCH)_e.o : - e_acosl.c e_atan2l.c e_exp2l.S e_expl.c e_fmodl.S e_log10l.S + e_acosl.c e_atan2l.c e_exp2l.S e_expl.S e_fmodl.S e_log10l.S e_log2l.S e_logl.S e_powl.S e_remainderl.S e_rem_pio2l.c e_scalbl.S e_sqrt.c e_sqrtf.c e_sqrtl.c ; diff --git a/src/system/libroot/posix/glibc/arch/x86_64/e_expl.S b/src/system/libroot/posix/glibc/arch/x86_64/e_expl.S new file mode 100644 index 0000000000..14dd29dcad --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/x86_64/e_expl.S @@ -0,0 +1,226 @@ +/* + * Written by J.T. Conklin . + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper . + */ + +/* + * The 8087 method for the exponential function is to calculate + * exp(x) = 2^(x log2(e)) + * after separating integer and fractional parts + * x log2(e) = i + f, |f| <= .5 + * 2^i is immediate but f needs to be precise for long double accuracy. + * Suppress range reduction error in computing f by the following. + * Separate x into integer and fractional parts + * x = xi + xf, |xf| <= .5 + * Separate log2(e) into the sum of an exact number c0 and small part c1. + * c0 + c1 = log2(e) to extra precision + * Then + * f = (c0 xi - i) + c0 xf + c1 x + * where c0 xi is exact and so also is (c0 xi - i). + * -- moshier@na-net.ornl.gov + */ + +#include + +#ifdef USE_AS_EXP10L +# define IEEE754_EXPL __ieee754_exp10l +# define EXPL_FINITE __exp10l_finite +# define FLDLOG fldl2t +#elif defined USE_AS_EXPM1L +# define IEEE754_EXPL __expm1l +# undef EXPL_FINITE +# define FLDLOG fldl2e +#else +# define IEEE754_EXPL __ieee754_expl +# define EXPL_FINITE __expl_finite +# define FLDLOG fldl2e +#endif + + .section .rodata.cst16,"aM",@progbits,16 + + .p2align 4 +#ifdef USE_AS_EXP10L + .type c0,@object +c0: .byte 0, 0, 0, 0, 0, 0, 0x9a, 0xd4, 0x00, 0x40 + .byte 0, 0, 0, 0, 0, 0 + ASM_SIZE_DIRECTIVE(c0) + .type c1,@object +c1: .byte 0x58, 0x92, 0xfc, 0x15, 0x37, 0x9a, 0x97, 0xf0, 0xef, 0x3f + .byte 0, 0, 0, 0, 0, 0 + ASM_SIZE_DIRECTIVE(c1) +#else + .type c0,@object +c0: .byte 0, 0, 0, 0, 0, 0, 0xaa, 0xb8, 0xff, 0x3f + .byte 0, 0, 0, 0, 0, 0 + ASM_SIZE_DIRECTIVE(c0) + .type c1,@object +c1: .byte 0x20, 0xfa, 0xee, 0xc2, 0x5f, 0x70, 0xa5, 0xec, 0xed, 0x3f + .byte 0, 0, 0, 0, 0, 0 + ASM_SIZE_DIRECTIVE(c1) +#endif +#ifndef USE_AS_EXPM1L + .type csat,@object +csat: .byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x0e, 0x40 + .byte 0, 0, 0, 0, 0, 0 + ASM_SIZE_DIRECTIVE(csat) + .type cmin,@object +cmin: .byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x1, 0 + .byte 0, 0, 0, 0, 0, 0 + ASM_SIZE_DIRECTIVE(cmin) +#endif + +#ifdef PIC +# define MO(op) op##(%rip) +#else +# define MO(op) op +#endif + + .text +ENTRY(IEEE754_EXPL) +#ifdef USE_AS_EXPM1L + movzwl 8+8(%rsp), %eax + xorb $0x80, %ah // invert sign bit (now 1 is "positive") + cmpl $0xc006, %eax // is num positive and exp >= 6 (number is >= 128.0)? + jae HIDDEN_JUMPTARGET (__expl) // (if num is denormal, it is at least >= 64.0) +#endif + fldt 8(%rsp) +/* I added the following ugly construct because expl(+-Inf) resulted + in NaN. The ugliness results from the bright minds at Intel. + For the i686 the code can be written better. + -- drepper@cygnus.com. */ + fxam /* Is NaN or +-Inf? */ +#ifdef USE_AS_EXPM1L + xorb $0x80, %ah + cmpl $0xc006, %eax + fstsw %ax + movb $0x45, %dh + jb 4f + + /* Below -64.0 (may be -NaN or -Inf). */ + andb %ah, %dh + cmpb $0x01, %dh + je 2f /* Is +-NaN, jump. */ + jmp 1f /* -large, possibly -Inf. */ + +4: /* In range -64.0 to 64.0 (may be +-0 but not NaN or +-Inf). */ + /* Test for +-0 as argument. */ + andb %ah, %dh + cmpb $0x40, %dh + je 2f + + /* Test for arguments that are small but not subnormal. */ + movzwl 8+8(%rsp), %eax + andl $0x7fff, %eax + cmpl $0x3fbf, %eax + jge 3f + /* Argument's exponent below -64; avoid spurious underflow if + normal. */ + cmpl $0x0001, %eax + jge 2f + /* Force underflow and return the argument, to avoid wrong signs + of zero results from the code below in some rounding modes. */ + fld %st + fmul %st + fstp %st + jmp 2f +#else + movzwl 8+8(%rsp), %eax + andl $0x7fff, %eax + cmpl $0x400d, %eax + jg 5f + cmpl $0x3fbc, %eax + jge 3f + /* Argument's exponent below -67, result rounds to 1. */ + fld1 + faddp + jmp 2f +5: /* Overflow, underflow or infinity or NaN as argument. */ + fstsw %ax + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + cmpb $0x01, %dh + je 2f /* Is +-NaN, jump. */ + /* Overflow or underflow; saturate. */ + fstp %st + fldt MO(csat) + andb $2, %ah + jz 3f + fchs +#endif +3: FLDLOG /* 1 log2(base) */ + fmul %st(1), %st /* 1 x log2(base) */ + /* Set round-to-nearest temporarily. */ + fstcw -4(%rsp) + movl $0xf3ff, %edx + andl -4(%rsp), %edx + movl %edx, -8(%rsp) + fldcw -8(%rsp) + frndint /* 1 i */ + fld %st(1) /* 2 x */ + frndint /* 2 xi */ + fldcw -4(%rsp) + fld %st(1) /* 3 i */ + fldt MO(c0) /* 4 c0 */ + fld %st(2) /* 5 xi */ + fmul %st(1), %st /* 5 c0 xi */ + fsubp %st, %st(2) /* 4 f = c0 xi - i */ + fld %st(4) /* 5 x */ + fsub %st(3), %st /* 5 xf = x - xi */ + fmulp %st, %st(1) /* 4 c0 xf */ + faddp %st, %st(1) /* 3 f = f + c0 xf */ + fldt MO(c1) /* 4 */ + fmul %st(4), %st /* 4 c1 * x */ + faddp %st, %st(1) /* 3 f = f + c1 * x */ + f2xm1 /* 3 2^(fract(x * log2(base))) - 1 */ +#ifdef USE_AS_EXPM1L + fstp %st(1) /* 2 */ + fscale /* 2 scale factor is st(1); base^x - 2^i */ + fxch /* 2 i */ + fld1 /* 3 1.0 */ + fscale /* 3 2^i */ + fld1 /* 4 1.0 */ + fsubrp %st, %st(1) /* 3 2^i - 1.0 */ + fstp %st(1) /* 2 */ + faddp %st, %st(1) /* 1 base^x - 1.0 */ +#else + fld1 /* 4 1.0 */ + faddp /* 3 2^(fract(x * log2(base))) */ + fstp %st(1) /* 2 */ + fscale /* 2 scale factor is st(1); base^x */ + fstp %st(1) /* 1 */ + /* Ensure underflow for tiny result. */ + fldt MO(cmin) /* 2 cmin */ + fld %st(1) /* 3 */ + fcomip %st(1), %st /* 2 */ + fstp %st /* 1 */ + jnc 6f + fld %st + fmul %st + fstp %st +#endif +6: fstp %st(1) /* 0 */ + jmp 2f +1: +#ifdef USE_AS_EXPM1L + /* For expm1l, only negative sign gets here. */ + fstp %st + fld1 + fchs +#else + testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ + fstp %st + fldz /* Set result to 0. */ +#endif +2: ret +END(IEEE754_EXPL) +#ifdef USE_AS_EXPM1L +libm_hidden_def (__expm1l) +weak_alias (__expm1l, expm1l) +#else +strong_alias (IEEE754_EXPL, EXPL_FINITE) +#endif diff --git a/src/system/libroot/posix/glibc/arch/x86_64/e_expl.c b/src/system/libroot/posix/glibc/arch/x86_64/e_expl.c deleted file mode 100644 index 1b299239e7..0000000000 --- a/src/system/libroot/posix/glibc/arch/x86_64/e_expl.c +++ /dev/null @@ -1,77 +0,0 @@ -/* - * Written by J.T. Conklin . - * Public domain. - * - * Adapted for `long double' by Ulrich Drepper . - */ - -/* - * The 8087 method for the exponential function is to calculate - * exp(x) = 2^(x log2(e)) - * after separating integer and fractional parts - * x log2(e) = i + f, |f| <= .5 - * 2^i is immediate but f needs to be precise for long double accuracy. - * Suppress range reduction error in computing f by the following. - * Separate x into integer and fractional parts - * x = xi + xf, |xf| <= .5 - * Separate log2(e) into the sum of an exact number c0 and small part c1. - * c0 + c1 = log2(e) to extra precision - * Then - * f = (c0 xi - i) + c0 xf + c1 x - * where c0 xi is exact and so also is (c0 xi - i). - * -- moshier@na-net.ornl.gov - */ - -#include "math_private.h" - -static const long double c0 = 1.44268798828125L; -static const long double c1 = 7.05260771340735992468e-6L; - -long double -__ieee754_expl (long double x) -{ - long double res; - -/* I added the following ugly construct because expl(+-Inf) resulted - in NaN. The ugliness results from the bright minds at Intel. - For the i686 the code can be written better. - -- drepper@cygnus.com. */ - asm ("fxam\n\t" /* Is NaN or +-Inf? */ - "fstsw %%ax\n\t" - "movb $0x45, %%dh\n\t" - "andb %%ah, %%dh\n\t" - "cmpb $0x05, %%dh\n\t" - "je 1f\n\t" /* Is +-Inf, jump. */ - "fldl2e\n\t" /* 1 log2(e) */ - "fmul %%st(1),%%st\n\t" /* 1 x log2(e) */ - "frndint\n\t" /* 1 i */ - "fld %%st(1)\n\t" /* 2 x */ - "frndint\n\t" /* 2 xi */ - "fld %%st(1)\n\t" /* 3 i */ - "fldt %2\n\t" /* 4 c0 */ - "fld %%st(2)\n\t" /* 5 xi */ - "fmul %%st(1),%%st\n\t" /* 5 c0 xi */ - "fsubp %%st,%%st(2)\n\t" /* 4 f = c0 xi - i */ - "fld %%st(4)\n\t" /* 5 x */ - "fsub %%st(3),%%st\n\t" /* 5 xf = x - xi */ - "fmulp %%st,%%st(1)\n\t" /* 4 c0 xf */ - "faddp %%st,%%st(1)\n\t" /* 3 f = f + c0 xf */ - "fldt %3\n\t" /* 4 */ - "fmul %%st(4),%%st\n\t" /* 4 c1 * x */ - "faddp %%st,%%st(1)\n\t" /* 3 f = f + c1 * x */ - "f2xm1\n\t" /* 3 2^(fract(x * log2(e))) - 1 */ - "fld1\n\t" /* 4 1.0 */ - "faddp\n\t" /* 3 2^(fract(x * log2(e))) */ - "fstp %%st(1)\n\t" /* 2 */ - "fscale\n\t" /* 2 scale factor is st(1); e^x */ - "fstp %%st(1)\n\t" /* 1 */ - "fstp %%st(1)\n\t" /* 0 */ - "jmp 2f\n\t" - "1:\ttestl $0x200, %%eax\n\t" /* Test sign. */ - "jz 2f\n\t" /* If positive, jump. */ - "fstp %%st\n\t" - "fldz\n\t" /* Set result to 0. */ - "2:\t\n" - : "=t" (res) : "0" (x), "m" (c0), "m" (c1) : "ax", "dx"); - return res; -} diff --git a/src/system/libroot/posix/glibc/arch/x86_64/s_expm1l.S b/src/system/libroot/posix/glibc/arch/x86_64/s_expm1l.S index eeffc3ebf5..af48f09082 100644 --- a/src/system/libroot/posix/glibc/arch/x86_64/s_expm1l.S +++ b/src/system/libroot/posix/glibc/arch/x86_64/s_expm1l.S @@ -1,89 +1,2 @@ -/* ix87 specific implementation of exp(x)-1. - Copyright (C) 1996,1997,2001,2002,2008,2009 Free Software Foundation, Inc. - This file is part of the GNU C Library. - Contributed by Ulrich Drepper , 1996. - Based on code by John C. Bowman . - Corrections by H.J. Lu (hjl@gnu.ai.mit.edu), 1997. - - 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, write to the Free - Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA - 02111-1307 USA. */ - - /* Using: e^x - 1 = 2^(x * log2(e)) - 1 */ - -#include - -#ifdef __ELF__ - .section .rodata -#else - .text -#endif - .align ALIGNARG(4) - ASM_TYPE_DIRECTIVE(minus1,@object) -minus1: .double -1.0 - ASM_SIZE_DIRECTIVE(minus1) - ASM_TYPE_DIRECTIVE(one,@object) -one: .double 1.0 - ASM_SIZE_DIRECTIVE(one) - ASM_TYPE_DIRECTIVE(l2e,@object) -l2e: .tfloat 1.442695040888963407359924681002 - ASM_SIZE_DIRECTIVE(l2e) - -#ifdef PIC -#define MO(op) op##(%rip) -#else -#define MO(op) op -#endif - - .text -ENTRY(__expm1l) - movzwl 8+8(%rsp), %eax // load sign bit and 15-bit exponent - xorb $0x80, %ah // invert sign bit (now 1 is "positive") - cmpl $0xc006, %eax // is num positive and exp >= 6 (number is >= 128.0)? - jae __expl@PLT // (if num is denormal, it is at least >= 64.0) - - fldt 8(%rsp) // x - fxam // Is NaN or +-Inf? - fstsw %ax - movb $0x45, %ch - andb %ah, %ch - cmpb $0x40, %ch - je 3f // If +-0, jump. - cmpb $0x05, %ch - je 2f // If +-Inf, jump. - - fldt MO(l2e) // log2(e) : x - fmulp // log2(e)*x - fld %st // log2(e)*x : log2(e)*x - frndint // int(log2(e)*x) : log2(e)*x - fsubr %st, %st(1) // int(log2(e)*x) : fract(log2(e)*x) - fxch // fract(log2(e)*x) : int(log2(e)*x) - f2xm1 // 2^fract(log2(e)*x)-1 : int(log2(e)*x) - fscale // 2^(log2(e)*x)-2^int(log2(e)*x) : int(log2(e)*x) - fxch // int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) - fldl MO(one) // 1 : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) - fscale // 2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) - fsubrl MO(one) // 1-2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) - fstp %st(1) // 1-2^int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) - fsubrp %st, %st(1) // 2^(log2(e)*x)-1 - ret - -2: testl $0x200, %eax // Test sign. - jz 3f // If positive, jump. - fstp %st - fldl MO(minus1) // Set result to -1.0. -3: ret -END(__expm1l) -libm_hidden_def (__expm1l) -weak_alias (__expm1l, expm1l) +#define USE_AS_EXPM1L +#include "e_expl.S"