libroot/glibc: Merge newer e_expl* files.
These versions (from ~2015 glibc) do not use some GCC-specific extensions that Clang's integrated assembler does not like.
This commit is contained in:
Augustin Cavalier
parent
b456e5d055
commit
834100c94f
@ -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
|
||||
;
|
||||
|
||||
226
src/system/libroot/posix/glibc/arch/x86_64/e_expl.S
Normal file
226
src/system/libroot/posix/glibc/arch/x86_64/e_expl.S
Normal file
@ -0,0 +1,226 @@
|
||||
/*
|
||||
* Written by J.T. Conklin <jtc@netbsd.org>.
|
||||
* Public domain.
|
||||
*
|
||||
* Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
|
||||
*/
|
||||
|
||||
/*
|
||||
* 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 <machine/asm.h>
|
||||
|
||||
#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
|
||||
@ -1,77 +0,0 @@
|
||||
/*
|
||||
* Written by J.T. Conklin <jtc@netbsd.org>.
|
||||
* Public domain.
|
||||
*
|
||||
* Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
|
||||
*/
|
||||
|
||||
/*
|
||||
* 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;
|
||||
}
|
||||
@ -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 <drepper@cygnus.com>, 1996.
|
||||
Based on code by John C. Bowman <bowman@ipp-garching.mpg.de>.
|
||||
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 <machine/asm.h>
|
||||
|
||||
#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"
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user