From patchwork Sun Sep 24 09:43:58 2017 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: "H.J. Lu" X-Patchwork-Id: 23120 Received: (qmail 28781 invoked by alias); 24 Sep 2017 09:44:05 -0000 Mailing-List: contact libc-alpha-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Unsubscribe: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: libc-alpha-owner@sourceware.org Delivered-To: mailing list libc-alpha@sourceware.org Received: (qmail 28764 invoked by uid 89); 24 Sep 2017 09:44:04 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-24.9 required=5.0 tests=AWL, BAYES_00, FREEMAIL_FROM, GIT_PATCH_0, GIT_PATCH_1, GIT_PATCH_2, GIT_PATCH_3, KAM_STOCKGEN, RCVD_IN_DNSWL_NONE, RCVD_IN_SORBS_SPAM, SPF_PASS autolearn=ham version=3.3.2 spammy= X-HELO: mail-oi0-f67.google.com X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20161025; h=x-gm-message-state:mime-version:in-reply-to:references:from:date :message-id:subject:to:cc; bh=oncA9HMxIkikA3CyAPBhlx19kr0jwO/UsLZml7zjQ+g=; b=qEEmm7P46D3pT/rWuGsyku7Cd9CVHUFilUiDMvCModWBCqIlKxalY+b5XbM4WNUWsQ nAiJG/bUIo6yw70ogJR6ahqhX9ro9QGcLW+fB2DzLlundqdBHtG2d89m9m2TwHl/2Imf rXvvjO8J28Rh3XRFYg8250XjiPICSkMAz5our3zYUTR42uV7SkZLa3qTRPRy7siLwjD+ Z1s2pFUXdCekBujojQOCEb5Gbt8j7nw4w4LYD2pimTRCoM9TYUmzxSrYc9/eBJqlyY4f jquP1ClKzUZyb6l2TO+ghuUUiliq16v0dA2jHIBOb3VkwFeGaoc4olHjDKQS1tXiGklU AfFg== X-Gm-Message-State: AHPjjUgmhR0dsl8cl0/FDjnndel/xlCkO9kVzx1JVS3OeV7giSs2J4iO 86dcmGqKNR7hEiCBGo760U6u5E5vqi3EoWweQ5A= X-Google-Smtp-Source: AOwi7QDgj8YbgzEhedyBu+tq7Ejl77Smir8YyqSaPhQW3RQH7r/GLj2xcCBk/JIMoKsZpJse7o1MRhK+GWxc0JlWsLQ= X-Received: by 10.202.59.212 with SMTP id i203mr2956504oia.103.1506246239327; Sun, 24 Sep 2017 02:43:59 -0700 (PDT) MIME-Version: 1.0 In-Reply-To: References: <59C1123F.9080003@arm.com> <59C23691.8080100@arm.com> From: "H.J. Lu" Date: Sun, 24 Sep 2017 02:43:58 -0700 Message-ID: Subject: Re: [PATCH 0/9] Optimized expf, exp2f, logf, log2f and powf To: Arjan van de Ven , Szabolcs Nagy Cc: GNU C Library , nd@arm.com On 9/23/17, H.J. Lu wrote: > > > On September 24, 2017 8:53:33 AM GMT+08:00, Arjan van de Ven > wrote: >>>>> >>>> >>>> now i added these to an nsz/math branch so others can test it. >>>> >>> >>> I created hjl/nsz/math branch to replace x86-64 assembly versions of >>> e_expf with generic >>> e_expf.c. I got >>> >> >> >> >>all results are <= 1ulp >> >>> >>> Test suite completed: >>> 356 test cases plus 352 tests for exception flags and >>> 352 tests for errno executed. >>> 6 errors occurred. >>> [hjl@gnu-efi-2 build-x86_64-linux]$ >>> >>> Arjan, have you measured performance assembly versions of expf vs >>generic expf? >>> >> >>the generic expf is faster than the ASM version; the generic version >>has an order 3 polyinomal >>while the ASM version has an order 4. >>In addition, the embedding of the errno stuff is a big deal (I already >>had that for exp() in my >>own project) since the exp[f]() code knows exactly when it's not needed >>at all (the fast path) >> >>so we really should move x86-64 over to this generic version if the >>1ulp is acceptable > > Agreed. I think 1 up is acceptable. I will prepare a patch > to do that. > Here is the patch on hjl/nsz/math branch with updated sysdeps/x86_64/fpu/libm-test-ulps. From f09e92b7a523a55dff2fdaf50699a78dcda63fde Mon Sep 17 00:00:00 2001 From: "H.J. Lu" Date: Sat, 23 Sep 2017 15:48:26 -0700 Subject: [PATCH] x86-64: Replace assembly versions of e_expf with generic e_expf.c This patch replaces x86-64 assembly versions of e_expf with generic e_expf.c. * sysdeps/x86_64/fpu/e_expf.S: Removed. * sysdeps/x86_64/fpu/multiarch/e_expf-fma.S: Likewise. * sysdeps/x86_64/fpu/w_expf.c: Likewise. * sysdeps/x86_64/fpu/libm-test-ulps: Updated for generic e_expf.c. * sysdeps/x86_64/fpu/multiarch/Makefile (CFLAGS-e_expf-fma.c): New. * sysdeps/x86_64/fpu/multiarch/e_expf-fma.c: New file. * sysdeps/x86_64/fpu/multiarch/e_expf.c (__redirect_ieee754_expf): Renamed to ... (__redirect_expf): This. (SYMBOL_NAME): Changed to expf. (__ieee754_expf): Renamed to ... (__expf): This. (__GI___expf): This. (__ieee754_expf): Add strong_alias. (__expf_finite): Likewise. (__expf): New. Include . --- sysdeps/x86_64/fpu/e_expf.S | 339 ------------------------------ sysdeps/x86_64/fpu/libm-test-ulps | 8 +- sysdeps/x86_64/fpu/multiarch/Makefile | 3 +- sysdeps/x86_64/fpu/multiarch/e_expf-fma.S | 182 ---------------- sysdeps/x86_64/fpu/multiarch/e_expf-fma.c | 3 + sysdeps/x86_64/fpu/multiarch/e_expf.c | 24 ++- sysdeps/x86_64/fpu/w_expf.c | 1 - 7 files changed, 30 insertions(+), 530 deletions(-) delete mode 100644 sysdeps/x86_64/fpu/e_expf.S delete mode 100644 sysdeps/x86_64/fpu/multiarch/e_expf-fma.S create mode 100644 sysdeps/x86_64/fpu/multiarch/e_expf-fma.c delete mode 100644 sysdeps/x86_64/fpu/w_expf.c diff --git a/sysdeps/x86_64/fpu/e_expf.S b/sysdeps/x86_64/fpu/e_expf.S deleted file mode 100644 index c3bf312c44..0000000000 --- a/sysdeps/x86_64/fpu/e_expf.S +++ /dev/null @@ -1,339 +0,0 @@ -/* Optimized __ieee754_expf function. - Copyright (C) 2012-2017 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, see - . */ - -#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 [0..K-1] - * 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) - /* Input: single precision x in %xmm0 */ - cvtss2sd %xmm0, %xmm1 /* Convert x to double precision */ - movd %xmm0, %ecx /* Copy x */ - movsd L(DP_KLN2)(%rip), %xmm2 /* DP K/log(2) */ - movsd L(DP_P2)(%rip), %xmm3 /* DP P2 */ - movl %ecx, %eax /* x */ - mulsd %xmm1, %xmm2 /* DP x*K/log(2) */ - andl $0x7fffffff, %ecx /* |x| */ - lea L(DP_T)(%rip), %rsi /* address of table T[j] */ - cmpl $0x42ad496b, %ecx /* |x|<125*log(2) ? */ - movsd L(DP_P3)(%rip), %xmm4 /* DP P3 */ - addsd L(DP_RS)(%rip), %xmm2 /* DP x*K/log(2)+RS */ - jae L(special_paths) - - /* Here if |x|<125*log(2) */ - cmpl $0x31800000, %ecx /* |x|<2^(-28) ? */ - jb L(small_arg) - - /* Main path: here if 2^(-28)<=|x|<125*log(2) */ - cvtsd2ss %xmm2, %xmm2 /* SP x*K/log(2)+RS */ - movd %xmm2, %eax /* bits of n*K+j with trash */ - subss L(SP_RS)(%rip), %xmm2 /* SP t=round(x*K/log(2)) */ - movl %eax, %edx /* n*K+j with trash */ - cvtss2sd %xmm2, %xmm2 /* DP t */ - andl $0x3f, %eax /* bits of j */ - mulsd L(DP_NLN2K)(%rip), %xmm2/* DP -t*log(2)/K */ - andl $0xffffffc0, %edx /* bits of n */ -#ifdef __AVX__ - vaddsd %xmm1, %xmm2, %xmm0 /* DP y=x-t*log(2)/K */ - vmulsd %xmm0, %xmm0, %xmm2 /* DP z=y*y */ -#else - addsd %xmm1, %xmm2 /* DP y=x-t*log(2)/K */ - movaps %xmm2, %xmm0 /* DP y */ - mulsd %xmm2, %xmm2 /* DP z=y*y */ -#endif - mulsd %xmm2, %xmm4 /* DP P3*z */ - addl $0x1fc0, %edx /* bits of n + SP exponent bias */ - mulsd %xmm2, %xmm3 /* DP P2*z */ - shll $17, %edx /* SP 2^n */ - addsd L(DP_P1)(%rip), %xmm4 /* DP P3*z+P1 */ - addsd L(DP_P0)(%rip), %xmm3 /* DP P2*z+P0 */ - movd %edx, %xmm1 /* SP 2^n */ - mulsd %xmm2, %xmm4 /* DP (P3*z+P1)*z */ - mulsd %xmm3, %xmm0 /* DP (P2*z+P0)*y */ - addsd %xmm4, %xmm0 /* DP P(y) */ - mulsd (%rsi,%rax,8), %xmm0 /* DP P(y)*T[j] */ - addsd (%rsi,%rax,8), %xmm0 /* DP T[j]*(P(y)+1) */ - cvtsd2ss %xmm0, %xmm0 /* SP T[j]*(P(y)+1) */ - mulss %xmm1, %xmm0 /* SP result=2^n*(T[j]*(P(y)+1)) */ - ret - - .p2align 4 -L(small_arg): - /* Here if 0<=|x|<2^(-28) */ - addss L(SP_ONE)(%rip), %xmm0 /* 1.0 + x */ - /* Return 1.0 with inexact raised, except for x==0 */ - ret - - .p2align 4 -L(special_paths): - /* 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 $0x7f800000, %ecx /* |x| is finite ? */ - jae L(arg_inf_or_nan) - - /* Here if |x|>under/overflow bound, and x is finite */ - testq %rax, %rax /* sign of x nonzero ? */ - je L(res_overflow) - - /* 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 - - .p2align 4 -L(arg_inf_or_nan): - /* Here if |x| is Inf or NAN */ - jne L(arg_nan) /* |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 - - .p2align 4 -L(arg_nan): - /* Here if |x| is NaN */ - addss %xmm0, %xmm0 /* Return x+x (raise invalid) */ - ret - - .p2align 4 -L(near_under_or_overflow): - /* Here if 125*log(2)<=|x|this bound, then result overflows */ - .long 0x42cff1b4 /* if x. */ - -#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 [0..K-1] - 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 - */ - - .section .text.fma,"ax",@progbits -ENTRY(__ieee754_expf_fma) - /* Input: single precision x in %xmm0 */ - vcvtss2sd %xmm0, %xmm0, %xmm1 /* Convert x to double precision */ - vmovd %xmm0, %ecx /* Copy x */ - vmovsd L(DP_KLN2)(%rip), %xmm2 /* DP K/log(2) */ - vfmadd213sd L(DP_RD)(%rip), %xmm1, %xmm2 /* DP x*K/log(2)+RD */ - vmovsd L(DP_P2)(%rip), %xmm3 /* DP P2 */ - movl %ecx, %eax /* x */ - andl $0x7fffffff, %ecx /* |x| */ - lea L(DP_T)(%rip), %rsi /* address of table T[j] */ - vmovsd L(DP_P3)(%rip), %xmm4 /* DP P3 */ - - cmpl $0x42ad496b, %ecx /* |x|<125*log(2) ? */ - jae L(special_paths_fma) - - /* Here if |x|<125*log(2) */ - cmpl $0x31800000, %ecx /* |x|<2^(-28) ? */ - jb L(small_arg_fma) - - /* Main path: here if 2^(-28)<=|x|<125*log(2) */ - /* %xmm2 = SP x*K/log(2)+RS */ - vmovd %xmm2, %eax - vsubsd L(DP_RD)(%rip), %xmm2, %xmm2 /* DP t=round(x*K/log(2)) */ - movl %eax, %edx /* n*K+j with trash */ - andl $0x3f, %eax /* bits of j */ - vmovsd (%rsi,%rax,8), %xmm5 /* T[j] */ - andl $0xffffffc0, %edx /* bits of n */ - - vfmadd132sd L(DP_NLN2K)(%rip), %xmm1, %xmm2 /* DP y=x-t*log(2)/K */ - vmulsd %xmm2, %xmm2, %xmm6 /* DP z=y*y */ - - - vfmadd213sd L(DP_P1)(%rip), %xmm6, %xmm4 /* DP P3*z + P1 */ - vfmadd213sd L(DP_P0)(%rip), %xmm6, %xmm3 /* DP P2*z+P0 */ - - addl $0x1fc0, %edx /* bits of n + SP exponent bias */ - shll $17, %edx /* SP 2^n */ - vmovd %edx, %xmm1 /* SP 2^n */ - - vmulsd %xmm6, %xmm4, %xmm4 /* DP (P3*z+P1)*z */ - - vfmadd213sd %xmm4, %xmm3, %xmm2 /* DP P(Y) (P2*z+P0)*y */ - vfmadd213sd %xmm5, %xmm5, %xmm2 /* DP T[j]*(P(y)+1) */ - vcvtsd2ss %xmm2, %xmm2, %xmm0 /* SP T[j]*(P(y)+1) */ - vmulss %xmm1, %xmm0, %xmm0 /* SP result=2^n*(T[j]*(P(y)+1)) */ - ret - - .p2align 4 -L(small_arg_fma): - /* Here if 0<=|x|<2^(-28) */ - vaddss L(SP_ONE)(%rip), %xmm0, %xmm0 /* 1.0 + x */ - /* Return 1.0 with inexact raised, except for x==0 */ - ret - - .p2align 4 -L(special_paths_fma): - /* 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 $0x7f800000, %ecx /* |x| is finite ? */ - jae L(arg_inf_or_nan_fma) - - /* Here if |x|>under/overflow bound, and x is finite */ - testl %eax, %eax /* sign of x nonzero ? */ - je L(res_overflow_fma) - - /* Here if -inf0) */ - vmovss L(SP_LARGE)(%rip), %xmm0/* load large value 2^100 */ - vmulss %xmm0, %xmm0, %xmm0 /* Return overflowed result (Inf or max normal) */ - ret - - .p2align 4 -L(arg_inf_or_nan_fma): - /* Here if |x| is Inf or NAN */ - jne L(arg_nan_fma) /* |x| is Inf ? */ - - /* Here if |x| is Inf */ - lea L(SP_INF_0)(%rip), %rdx /* depending on sign of x: */ - vmovss (%rdx,%rax,4), %xmm0 /* return zero or Inf */ - ret - - .p2align 4 -L(arg_nan_fma): - /* Here if |x| is NaN */ - vaddss %xmm0, %xmm0, %xmm0 /* Return x+x (raise invalid) */ - ret - - .p2align 4 -L(near_under_or_overflow_fma): - /* Here if 125*log(2)<=|x| diff --git a/sysdeps/x86_64/fpu/multiarch/e_expf-fma.c b/sysdeps/x86_64/fpu/multiarch/e_expf-fma.c new file mode 100644 index 0000000000..4e01cd6a82 --- /dev/null +++ b/sysdeps/x86_64/fpu/multiarch/e_expf-fma.c @@ -0,0 +1,3 @@ +#define __expf __expf_fma + +#include diff --git a/sysdeps/x86_64/fpu/multiarch/e_expf.c b/sysdeps/x86_64/fpu/multiarch/e_expf.c index 096209857e..3cdbf6cf2a 100644 --- a/sysdeps/x86_64/fpu/multiarch/e_expf.c +++ b/sysdeps/x86_64/fpu/multiarch/e_expf.c @@ -16,11 +16,25 @@ License along with the GNU C Library; if not, see . */ -extern float __redirect_ieee754_expf (float); +extern float __redirect_expf (float); -#define SYMBOL_NAME ieee754_expf +#define SYMBOL_NAME expf #include "ifunc-fma.h" -libc_ifunc_redirected (__redirect_ieee754_expf, __ieee754_expf, - IFUNC_SELECTOR ()); -strong_alias (__ieee754_expf, __expf_finite) +libc_ifunc_redirected (__redirect_expf, __expf, IFUNC_SELECTOR ()); + +#ifdef SHARED +__hidden_ver1 (__expf, __GI___expf, __redirect_expf) + __attribute__ ((visibility ("hidden"))); + +# include +versioned_symbol (libm, __expf, expf, GLIBC_2_27); +#else +weak_alias (__expf, expf) +#endif + +strong_alias (__expf, __ieee754_expf) +strong_alias (__expf, __expf_finite) + +#define __expf __expf_sse2 +#include diff --git a/sysdeps/x86_64/fpu/w_expf.c b/sysdeps/x86_64/fpu/w_expf.c deleted file mode 100644 index b5fe164520..0000000000 --- a/sysdeps/x86_64/fpu/w_expf.c +++ /dev/null @@ -1 +0,0 @@ -#include -- 2.13.5