diff --git a/sysdeps/x86_64/fpu/e_expf.S b/sysdeps/x86_64/fpu/e_expf.S new file mode 100755 index 0000000..c0d0e18 --- /dev/null +++ b/sysdeps/x86_64/fpu/e_expf.S @@ -0,0 +1,345 @@ +/* Optimized __ieee754_expf function. + Copyright (C) 2012 Free Software Foundation, Inc. + Contributed by Intel Corporation. + This file is part of the GNU C Library. + + 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. */ + +#include + +/* Short algorithm description: + * + * Let K = 64 (table size). + * e^x = 2^(x/log(2)) = 2^n * T[j] * (1 + P(y)) + * where + * x = m*log(2)/K + y, y in [0.0..log(2)/K] + * m = n*K + j, m,n,j - signed integer, j in [-K/2..K/2] + * values of 2^(j/K) are tabulated as T[j]. + * + * P(y) is a minimax polynomial approximation of expf(x)-1 + * on small interval [0.0..log(2)/K]. + * + * P(y) = P3*y*y*y*y + P2*y*y*y + P1*y*y + P0*y, calculated as + * z = y*y; P(y) = (P3*z + P1)*z + (P2*z + P0)*y + * + * Special cases: + * expf(NaN) = NaN + * expf(+INF) = +INF + * expf(-INF) = 0 + * expf(x) = 1 for subnormals + * for finite argument, only expf(0)=1 is exact + * expf(x) overflows if x>88.7228317260742190 + * expf(x) underflows if x<-103.972076416015620 + */ + + .text +ENTRY(__ieee754_expf) + movd %xmm0, %eax /* Copy argument x */ + movss %xmm0, -8(%rsp) /* Save argument in current frame */ + movl %eax, %ecx + andl $2147483647, %ecx /* Clear sign bit */ + cmpl $1118652779, %ecx /* |x|<125*log(2) ? */ + jae L(label_1_6) +L(label_1_2): + /* Here if |x|<125*log(2) */ + cmpl $830472192, %ecx /* |x|<2^(-28) ? */ + jae L(label_1_4) +L(label_1_3): + /* Here if 0<=|x|<2^(-28) */ + movss L(SP_ONE)(%rip), %xmm1 /* 1.0 */ + addss %xmm1, %xmm0 /* 1.0 + x */ + /* Return 1.0 with inexact raised, except for x==0 */ + ret + +L(label_1_4): + /* Here if 2^(-28)<=|x|<125*log(2) */ + cvtss2sd -8(%rsp), %xmm3 /* Load x converted to double precision */ + movaps %xmm3, %xmm0 /* DP x */ + lea L(DP_T)(%rip), %rcx /* address of table T[j] */ + movsd L(DP_P2)(%rip), %xmm2 /* DP P2 */ + mulsd L(DP_KLN2)(%rip), %xmm0 /* DP x*K/log(2) */ + addsd L(DP_RS)(%rip), %xmm0 /* DP x*K/log(2)+RS */ + cvtsd2ss %xmm0, %xmm0 /* SP x*K/log(2)+RS */ + movss %xmm0, -16(%rsp) /* save SP x*K/log(2)+RS */ + movss -16(%rsp), %xmm1 /* load SP x*K/log(2)+RS */ + movsd L(DP_P3)(%rip), %xmm0 /* DP P3 */ + movl -16(%rsp), %edx /* n*K+j */ + movl %edx, %eax /* n*K+j */ + shll $26, %eax /* bits of j */ + subss L(SP_RS)(%rip), %xmm1 /* SP t=round(x*K/log(2)) */ + cvtss2sd %xmm1, %xmm1 /* DP t */ + mulsd L(DP_NLN2K)(%rip), %xmm1 /* DP -t*log(2)/K */ + sarl $26, %eax /* j */ + addsd %xmm1, %xmm3 /* DP y=x-t*log(2)/K */ + movaps %xmm3, %xmm4 /* DP y */ + mulsd %xmm3, %xmm4 /* DP z=y*y */ + mulsd %xmm4, %xmm0 /* DP P3*z */ + movslq %eax, %rax /* sign extend j (for table lookup index) */ + addsd L(DP_P1)(%rip), %xmm0 /* DP P3*z+P1 */ + mulsd %xmm4, %xmm0 /* DP (P3*z+P1)*z */ + subl %eax, %edx /* n*K */ + mulsd %xmm2, %xmm4 /* DP P2*z */ + shrl $6, %edx /* n */ + addsd L(DP_P0)(%rip), %xmm4 /* DP P2*z+P0 */ + mulsd %xmm3, %xmm4 /* DP (P2*z+P0)*y */ + shll $23, %edx /* n shifted to SP exponent field */ + addsd %xmm4, %xmm0 /* DP P(y) */ + mulsd 256(%rcx,%rax,8),%xmm0 /* DP P(y)*T[j] */ + addl $1065353216, %edx /* SP 2^n (added SP bias) */ + movl %edx, -24(%rsp) /* save SP 2^n */ + cvtss2sd -24(%rsp), %xmm5 /* load DP 2^n */ + addsd 256(%rcx,%rax,8),%xmm0 /* DP T[j]*(P(y)+1) */ + mulsd %xmm5, %xmm0 /* DP result=2^n*(T[j]*(P(y)+1)) */ + cvtsd2ss %xmm0, %xmm0 /* convert result to single precision */ + ret + +L(label_1_6): + /* Here if 125*log(2)<=|x| */ + shrl $31, %eax /* Get sign bit of x, and depending on it: */ + lea L(SP_RANGE)(%rip), %rdx /* load over/underflow bound */ + cmpl (%rdx, %rax, 4), %ecx /* |x|under/overflow bound */ + cmpl $2139095040, %ecx /* |x| is finite ? */ + jae L(label_1_12) + + /* Here if |x|>under/overflow bound, and x is finite */ + testq %rax, %rax /* sign of x nonzero ? */ + je L(label_1_10) + + /* Here if -inf0) */ + movss L(SP_LARGE)(%rip), %xmm0 /* load large value 2^100 */ + mulss %xmm0, %xmm0 /* Return overflowed result (Inf or max normal) */ + ret + +L(label_1_12): + /* Here if |x| is Inf or NAN */ + jne L(label_1_14) /* |x| is Inf ? */ + + /* Here if |x| is Inf */ + lea L(SP_INF_0)(%rip), %rdx /* depending on sign of x: */ + movss (%rdx, %rax, 4), %xmm0 /* return zero or Inf */ + ret + +L(label_1_14): + /* Here if |x| is NaN */ + movss -8(%rsp), %xmm0 /* Load x */ + addss %xmm0, %xmm0 /* Return x+x (raise invalid) */ + ret + +L(label_1_16): + /* Here if 125*log(2)<=|x|this bound, then result overflows */ + .long 1120924084 /* if x