From patchwork Thu Apr 9 19:59:42 2020 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Adhemerval Zanella Netto X-Patchwork-Id: 38793 Return-Path: X-Original-To: libc-alpha@sourceware.org Delivered-To: libc-alpha@sourceware.org Received: from mail-qk1-x730.google.com (mail-qk1-x730.google.com [IPv6:2607:f8b0:4864:20::730]) by sourceware.org (Postfix) with ESMTPS id D1C9B385B835 for ; Thu, 9 Apr 2020 19:59:52 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.3.2 sourceware.org D1C9B385B835 Received: by mail-qk1-x730.google.com with SMTP id 130so18053qke.4 for ; Thu, 09 Apr 2020 12:59:52 -0700 (PDT) X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20161025; h=x-gm-message-state:from:to:cc:subject:date:message-id:in-reply-to :references; bh=Dss1JJURnb+yMjHYsDB+KmrTIRmQls161ECw5TT6lSA=; b=qRmMDu7dltxVkCOCaME6iv4l9MVWWcoO0+aFWB7CxpwlBnKpNdtkxmYTAeflQ5fmSH 483kldRSfHPz1mAxNl5gge2R/16iIaQF3ArhFs8ZD6eLw9HqUrOKIWZTI10+3OrQGR1L kbCRE/u01g3puvu1zqkRjt+ixIf0UXglz/LROh5NplJIxyGW719xd+Y53G9wtwYAgE/V vaSQhAzETKH7QDAngEAPzEz1vMd+HhM/1URg8bBCtiDWv0lZflrFK8I9LoMgQ4nlfRHX lDr4S+f9hveBlwBfBwp6b4niwzig6DA1sPZxfHZ8Y60Q9U3htyTGDFmMlGQCBPe04M5Q mh1w== X-Gm-Message-State: AGi0PuZv1axRrb/9kbVD5rmcUEdO9RCga2UfcvE7bjtSDRGFDFceJsoV 40lqSqhXE92ZxM4tW0ntLGSIaDZfnnx1pw== X-Google-Smtp-Source: APiQypIal8hLmjyIZQVy9WHgDZPG9s2eYP5mx/UoUW00cGIEpeMY3z4IiyN8OWmDWmvnVEWZTmKtng== X-Received: by 2002:a37:b17:: with SMTP id 23mr582511qkl.136.1586462391804; Thu, 09 Apr 2020 12:59:51 -0700 (PDT) Received: from localhost.localdomain ([177.194.48.209]) by smtp.googlemail.com with ESMTPSA id f7sm8605434qkh.25.2020.04.09.12.59.50 (version=TLS1_3 cipher=TLS_AES_256_GCM_SHA384 bits=256/256); Thu, 09 Apr 2020 12:59:51 -0700 (PDT) From: Adhemerval Zanella To: libc-alpha@sourceware.org Subject: [PATCH 2/5] math: Optimized generic exp10f with wrappers Date: Thu, 9 Apr 2020 16:59:42 -0300 Message-Id: <20200409195945.10393-2-adhemerval.zanella@linaro.org> X-Mailer: git-send-email 2.17.1 In-Reply-To: <20200409195945.10393-1-adhemerval.zanella@linaro.org> References: <20200409195945.10393-1-adhemerval.zanella@linaro.org> X-Spam-Status: No, score=-25.6 required=5.0 tests=BAYES_00, DKIM_SIGNED, DKIM_VALID, DKIM_VALID_AU, DKIM_VALID_EF, GIT_PATCH_0, GIT_PATCH_1, GIT_PATCH_2, GIT_PATCH_3, KAM_SHORT, RCVD_IN_DNSWL_NONE, SPF_HELO_NONE, SPF_PASS, TXREP autolearn=ham autolearn_force=no version=3.4.2 X-Spam-Checker-Version: SpamAssassin 3.4.2 (2018-09-13) on server2.sourceware.org X-BeenThere: libc-alpha@sourceware.org X-Mailman-Version: 2.1.29 Precedence: list List-Id: Libc-alpha mailing list List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Thu, 09 Apr 2020 19:59:55 -0000 From: Paul Zimmermann It is inspired by expf and reuses its tables and internal functions. The error checks are inlined and errno setting is in separate tail called functions, but the wrappers are kept in this patch to handle the _LIB_VERSION==_SVID_ case. Double precision arithmetics is used which is expected to be faster on most targets (including soft-float) than using single precision and it is easier to get good precision result with it. Result for x86_64 (i7-4790K CPU @ 4.00GHz) are: Before new code: "exp10f": { "": { "duration": 1.00923e+09, "iterations": 5.4832e+07, "max": 249.361, "min": 12.5, "mean": 18.4059 } With new code: "exp10f": { "": { "duration": 9.7639e+08, "iterations": 8.1056e+07, "max": 563.323, "min": 11.54, "mean": 12.0459 } Result for aarch64 (A72 @ 2GHz) are: Before new code: "exp10f": { "": { "duration": 1.00923e+09, "iterations": 5.4832e+07, "max": 249.361, "min": 12.5, "mean": 18.4059 } With new code: "exp10f": { "": { "duration": 9.7639e+08, "iterations": 8.1056e+07, "max": 563.323, "min": 11.54, "mean": 12.0459 } Checked on x86_64-linux-gnu, powerpc64le-linux-gnu, aarch64-linux-gnu, and sparc64-linux-gnu. --- math/e_exp10f.c | 32 ----- sysdeps/ieee754/flt-32/e_exp10f.c | 198 +++++++++++++++++++++++++++ sysdeps/ieee754/flt-32/math_config.h | 2 +- 3 files changed, 199 insertions(+), 33 deletions(-) delete mode 100644 math/e_exp10f.c create mode 100644 sysdeps/ieee754/flt-32/e_exp10f.c diff --git a/math/e_exp10f.c b/math/e_exp10f.c deleted file mode 100644 index 93c41d00a6..0000000000 --- a/math/e_exp10f.c +++ /dev/null @@ -1,32 +0,0 @@ -/* Copyright (C) 1998-2020 Free Software Foundation, Inc. - This file is part of the GNU C Library. - Contributed by Ulrich Drepper , 1998. - - 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 -#include -#include - -float -__ieee754_exp10f (float arg) -{ - /* The argument to exp needs to be calculated in enough precision - that the fractional part has as much precision as float, in - addition to the bits in the integer part; using double ensures - this. */ - return __ieee754_exp (M_LN10 * arg); -} -libm_alias_finite (__ieee754_exp10f, __exp10f) diff --git a/sysdeps/ieee754/flt-32/e_exp10f.c b/sysdeps/ieee754/flt-32/e_exp10f.c new file mode 100644 index 0000000000..034a9e364a --- /dev/null +++ b/sysdeps/ieee754/flt-32/e_exp10f.c @@ -0,0 +1,198 @@ +/* Single-precision 10^x function. + Copyright (C) 2020 Free Software Foundation, Inc. + 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 +#include +#include +#include +#include +#include "math_config.h" + +/* + EXP2F_TABLE_BITS 5 + EXP2F_POLY_ORDER 3 + + Max. ULP error: 0.502 (normal range, nearest rounding). + Max. relative error: 2^-33.240 (before rounding to float). + Wrong count: 169839. + Non-nearest ULP error: 1 (rounded ULP error). + + Detailed error analysis (for EXP2F_TABLE_BITS=3 thus N=32): + + - We first compute z = RN(InvLn10N * x) where + + InvLn10N = N*log(10)/log(2) * (1 + theta1) with |theta1| < 2^-54.150 + + since z is rounded to nearest: + + z = InvLn10N * x * (1 + theta2) with |theta2| < 2^-53 + + thus z = N*log(10)/log(2) * x * (1 + theta3) with |theta3| < 2^-52.463 + + - Since |x| < 38 in the main branch, we deduce: + + z = N*log(10)/log(2) * x + theta4 with |theta4| < 2^-40.483 + + - We then write z = k + r where k is an integer and |r| <= 0.5 (exact) + + - We thus have + + x = z*log(2)/(N*log(10)) - theta4*log(2)/(N*log(10)) + = z*log(2)/(N*log(10)) + theta5 with |theta5| < 2^-47.215 + + 10^x = 2^(k/N) * 2^(r/N) * 10^theta5 + = 2^(k/N) * 2^(r/N) * (1 + theta6) with |theta6| < 2^-46.011 + + - Then 2^(k/N) is approximated by table lookup, the maximal relative error + is for (k%N) = 5, with + + s = 2^(i/N) * (1 + theta7) with |theta7| < 2^-53.249 + + - 2^(r/N) is approximated by a degree-3 polynomial, where the maximal + mathematical relative error is 2^-33.243. + + - For C[0] * r + C[1], assuming no FMA is used, since |r| <= 0.5 and + |C[0]| < 1.694e-6, |C[0] * r| < 8.47e-7, and the absolute error on + C[0] * r is bounded by 1/2*ulp(8.47e-7) = 2^-74. Then we add C[1] with + |C[1]| < 0.000235, thus the absolute error on C[0] * r + C[1] is bounded + by 2^-65.994 (z is bounded by 0.000236). + + - For r2 = r * r, since |r| <= 0.5, we have |r2| <= 0.25 and the absolute + error is bounded by 1/4*ulp(0.25) = 2^-56 (the factor 1/4 is because |r2| + cannot exceed 1/4, and on the left of 1/4 the distance between two + consecutive numbers is 1/4*ulp(1/4)). + + - For y = C[2] * r + 1, assuming no FMA is used, since |r| <= 0.5 and + |C[2]| < 0.0217, the absolute error on C[2] * r is bounded by 2^-60, + and thus the absolute error on C[2] * r + 1 is bounded by 1/2*ulp(1)+2^60 + < 2^-52.988, and |y| < 1.01085 (the error bound is better if a FMA is + used). + + - for z * r2 + y: the absolute error on z is bounded by 2^-65.994, with + |z| < 0.000236, and the absolute error on r2 is bounded by 2^-56, with + r2 < 0.25, thus |z*r2| < 0.000059, and the absolute error on z * r2 + (including the rounding error) is bounded by: + + 2^-65.994 * 0.25 + 0.000236 * 2^-56 + 1/2*ulp(0.000059) < 2^-66.429. + + Now we add y, with |y| < 1.01085, and error on y bounded by 2^-52.988, + thus the absolute error is bounded by: + + 2^-66.429 + 2^-52.988 + 1/2*ulp(1.01085) < 2^-51.993 + + - Now we convert the error on y into relative error. Recall that y + approximates 2^(r/N), for |r| <= 0.5 and N=32. Thus 2^(-0.5/32) <= y, + and the relative error on y is bounded by + + 2^-51.993/2^(-0.5/32) < 2^-51.977 + + - Taking into account the mathematical relative error of 2^-33.243 we have: + + y = 2^(r/N) * (1 + theta8) with |theta8| < 2^-33.242 + + - Since we had s = 2^(k/N) * (1 + theta7) with |theta7| < 2^-53.249, + after y = y * s we get y = 2^(k/N) * 2^(r/N) * (1 + theta9) + with |theta9| < 2^-33.241 + + - Finally, taking into account the error theta6 due to the rounding error on + InvLn10N, and when multiplying InvLn10N * x, we get: + + y = 10^x * (1 + theta10) with |theta10| < 2^-33.240 + + - Converting into binary64 ulps, since y < 2^53*ulp(y), the error is at most + 2^19.76 ulp(y) + + - If the result is a binary32 value in the normal range (i.e., >= 2^-126), + then the error is at most 2^-9.24 ulps, i.e., less than 0.00166 (in the + subnormal range, the error in ulps might be larger). + + Note that this bound is tight, since for x=-0x25.54ac0p0 the final value of + y (before conversion to float) differs from 879853 ulps from the correctly + rounded value, and 879853 ~ 2^19.7469. */ + +#define N (1 << EXP2F_TABLE_BITS) + +#define InvLn10N (0x3.5269e12f346e2p0 * N) /* log(10)/log(2) to nearest */ +#define T __exp2f_data.tab +#define C __exp2f_data.poly_scaled +#define SHIFT __exp2f_data.shift + +static inline uint32_t +top13 (float x) +{ + return asuint (x) >> 19; +} + +float +__ieee754_exp10f (float x) +{ + uint32_t abstop; + uint64_t ki, t; + double kd, xd, z, r, r2, y, s; + + xd = (double) x; + abstop = top13 (x) & 0xfff; /* Ignore sign. */ + if (__glibc_unlikely (abstop >= top13 (38.0f))) + { + /* |x| >= 38 or x is nan. */ + if (asuint (x) == asuint (-INFINITY)) + return 0.0f; + if (abstop >= top13 (INFINITY)) + return x + x; + /* 0x26.8826ap0 is the largest value such that 10^x < 2^128. */ + if (x > 0x26.8826ap0f) + return __math_oflowf (0); + /* -0x2d.278d4p0 is the smallest value such that 10^x > 2^-150. */ + if (x < -0x2d.278d4p0f) + return __math_uflowf (0); +#if WANT_ERRNO_UFLOW + if (x < -0x2c.da7cfp0) + return __math_may_uflowf (0); +#endif + /* the smallest value such that 10^x >= 2^-126 (normal range) + is x = -0x25.ee060p0 */ + /* we go through here for 2014929 values out of 2060451840 + (not counting NaN and infinities, i.e., about 0.1% */ + } + + /* x*N*Ln10/Ln2 = k + r with r in [-1/2, 1/2] and int k. */ + z = InvLn10N * xd; + /* |xd| < 38 thus |z| < 1216 */ +#if TOINT_INTRINSICS + kd = roundtoint (z); + ki = converttoint (z); +#else +# define SHIFT __exp2f_data.shift + kd = math_narrow_eval ((double) (z + SHIFT)); /* Needs to be double. */ + ki = asuint64 (kd); + kd -= SHIFT; +#endif + r = z - kd; + + /* 10^x = 10^(k/N) * 10^(r/N) ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ + t = T[ki % N]; + t += ki << (52 - EXP2F_TABLE_BITS); + s = asdouble (t); + z = C[0] * r + C[1]; + r2 = r * r; + y = C[2] * r + 1; + y = z * r2 + y; + y = y * s; + return (float) y; +} +libm_alias_finite (__ieee754_exp10f, __exp10f) diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h index bf79274ce7..4817e500e1 100644 --- a/sysdeps/ieee754/flt-32/math_config.h +++ b/sysdeps/ieee754/flt-32/math_config.h @@ -109,7 +109,7 @@ attribute_hidden float __math_may_uflowf (uint32_t); attribute_hidden float __math_divzerof (uint32_t); attribute_hidden float __math_invalidf (float); -/* Shared between expf, exp2f and powf. */ +/* Shared between expf, exp2f, exp10f, and powf. */ #define EXP2F_TABLE_BITS 5 #define EXP2F_POLY_ORDER 3 extern const struct exp2f_data