From e9e23e58ca7404802d4f1e752cef38861771df2d Mon Sep 17 00:00:00 2001 From: Pierre Blanchard Date: Mon, 11 Apr 2022 14:42:01 +0100 Subject: pl/math: Add scalar log10f - Scalar log10f simply implemented as log10(x):=log(x)*invlog10. - The implementation of the log is similar to that of math/. - The maximum measured ULP error is about 0.80ulp. --- pl/math/log10f.c | 90 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 90 insertions(+) create mode 100644 pl/math/log10f.c (limited to 'pl/math/log10f.c') diff --git a/pl/math/log10f.c b/pl/math/log10f.c new file mode 100644 index 0000000..79f5d12 --- /dev/null +++ b/pl/math/log10f.c @@ -0,0 +1,90 @@ +/* + * Single-precision log10 function. + * + * Copyright (c) 2022, Arm Limited. + * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + */ + +#include "math_config.h" +#include +#include + +/* Data associated to logf: + + LOGF_TABLE_BITS = 4 + LOGF_POLY_ORDER = 4 + + ULP error: 0.818 (nearest rounding.) + Relative error: 1.957 * 2^-26 (before rounding.). */ + +#define T __logf_data.tab +#define A __logf_data.poly +#define Ln2 __logf_data.ln2 +#define InvLn10 __logf_data.invln10 +#define N (1 << LOGF_TABLE_BITS) +#define OFF 0x3f330000 + +/* This naive implementation of log10f mimics that of log + then simply scales the result by 1/log(10) to switch from base e to + base 10. Hence, most computations are carried out in double precision. + Scaling before rounding to single precision is both faster and more accurate. + + ULP error: 0.797 ulp (nearest rounding.). */ +float +log10f (float x) +{ + /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */ + double_t z, r, r2, y, y0, invc, logc; + uint32_t ix, iz, tmp; + int k, i; + + ix = asuint (x); +#if WANT_ROUNDING + /* Fix sign of zero with downward rounding when x==1. */ + if (unlikely (ix == 0x3f800000)) + return 0; +#endif + if (unlikely (ix - 0x00800000 >= 0x7f800000 - 0x00800000)) + { + /* x < 0x1p-126 or inf or nan. */ + if (ix * 2 == 0) + return __math_divzerof (1); + if (ix == 0x7f800000) /* log(inf) == inf. */ + return x; + if ((ix & 0x80000000) || ix * 2 >= 0xff000000) + return __math_invalidf (x); + /* x is subnormal, normalize it. */ + ix = asuint (x * 0x1p23f); + ix -= 23 << 23; + } + + /* x = 2^k z; where z is in range [OFF,2*OFF] and exact. + The range is split into N subintervals. + The ith subinterval contains z and c is near its center. */ + tmp = ix - OFF; + i = (tmp >> (23 - LOGF_TABLE_BITS)) % N; + k = (int32_t) tmp >> 23; /* arithmetic shift. */ + iz = ix - (tmp & 0x1ff << 23); + invc = T[i].invc; + logc = T[i].logc; + z = (double_t) asfloat (iz); + + /* log(x) = log1p(z/c-1) + log(c) + k*Ln2. */ + r = z * invc - 1; + y0 = logc + (double_t) k * Ln2; + + /* Pipelined polynomial evaluation to approximate log1p(r). */ + r2 = r * r; + y = A[1] * r + A[2]; + y = A[0] * r2 + y; + y = y * r2 + (y0 + r); + + /* Multiply by 1/log(10). */ + y = y * InvLn10; + + return eval_as_float (y); +} +#if USE_GLIBC_ABI +strong_alias (log10f, __log10f_finite) +hidden_alias (log10f, __ieee754_log10f) +#endif -- cgit v1.2.3 From 1bca1a541cce13c352296acd5dfa16160fc27bc9 Mon Sep 17 00:00:00 2001 From: Joe Ramsay Date: Thu, 15 Dec 2022 13:27:31 +0000 Subject: pl/math: Auto-generate mathbench and ulp headers Instead of maintaining three separate lists of routines, which are cumbersome and prone to merge conflicts, we provide a new macro, PL_SIG, which by some preprocessor machinery outputs the lists in the required format (macro formats have been changed very slightly to make the generation simpler). Only routines with simple signatures are handled - binary functions still need mathbench wrappers defined manually. As well, routines with non-standard references (i.e. powi/powk) still need entries and wrappers manually defined. --- pl/math/log10f.c | 6 ++++++ 1 file changed, 6 insertions(+) (limited to 'pl/math/log10f.c') diff --git a/pl/math/log10f.c b/pl/math/log10f.c index 79f5d12..ea67b4b 100644 --- a/pl/math/log10f.c +++ b/pl/math/log10f.c @@ -6,6 +6,7 @@ */ #include "math_config.h" +#include "pl_sig.h" #include #include @@ -84,7 +85,12 @@ log10f (float x) return eval_as_float (y); } + +// clang-format off #if USE_GLIBC_ABI strong_alias (log10f, __log10f_finite) hidden_alias (log10f, __ieee754_log10f) #endif + +PL_SIG (S, F, 1, log10, 0.01, 11.1) + // clang-format on -- cgit v1.2.3 From ecb1c6f6ea7872645cb4c26514d5f64815b61a1b Mon Sep 17 00:00:00 2001 From: Joe Ramsay Date: Thu, 15 Dec 2022 13:27:39 +0000 Subject: pl/math: Move ULP limits to routine source files Introduces a new set of macros and Make rules for mechanically generating a list of ULP limits for each routine, to be consumed by runulp.sh. This removes the need to maintain long lists of thresholds in runulp.sh. --- pl/math/log10f.c | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) (limited to 'pl/math/log10f.c') diff --git a/pl/math/log10f.c b/pl/math/log10f.c index ea67b4b..84db420 100644 --- a/pl/math/log10f.c +++ b/pl/math/log10f.c @@ -5,11 +5,13 @@ * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception */ -#include "math_config.h" -#include "pl_sig.h" #include #include +#include "math_config.h" +#include "pl_sig.h" +#include "pl_test.h" + /* Data associated to logf: LOGF_TABLE_BITS = 4 @@ -86,11 +88,5 @@ log10f (float x) return eval_as_float (y); } -// clang-format off -#if USE_GLIBC_ABI -strong_alias (log10f, __log10f_finite) -hidden_alias (log10f, __ieee754_log10f) -#endif - PL_SIG (S, F, 1, log10, 0.01, 11.1) - // clang-format on +PL_TEST_ULP (log10f, 0.30) -- cgit v1.2.3 From 202e46317ee8983516b6413066a57bd624ffa044 Mon Sep 17 00:00:00 2001 From: Joe Ramsay Date: Thu, 15 Dec 2022 13:28:06 +0000 Subject: pl/math: Move test intervals to routine source files To conclude the work on simplifying the runulp.sh script, a new macro has been introduced to specify the intervals in which a routine should be tested in the routine source. This is eventually consumed by runulp.sh. --- pl/math/log10f.c | 5 +++++ 1 file changed, 5 insertions(+) (limited to 'pl/math/log10f.c') diff --git a/pl/math/log10f.c b/pl/math/log10f.c index 84db420..32de42f 100644 --- a/pl/math/log10f.c +++ b/pl/math/log10f.c @@ -90,3 +90,8 @@ log10f (float x) PL_SIG (S, F, 1, log10, 0.01, 11.1) PL_TEST_ULP (log10f, 0.30) +PL_TEST_INTERVAL (log10f, 0, 0xffff0000, 10000) +PL_TEST_INTERVAL (log10f, 0x1p-127, 0x1p-26, 50000) +PL_TEST_INTERVAL (log10f, 0x1p-26, 0x1p3, 50000) +PL_TEST_INTERVAL (log10f, 0x1p-4, 0x1p4, 50000) +PL_TEST_INTERVAL (log10f, 0, inf, 50000) -- cgit v1.2.3 From 0f87f607b976820ef41fe64d004fe67dc7af8236 Mon Sep 17 00:00:00 2001 From: Joe Ramsay Date: Thu, 5 Jan 2023 11:56:20 +0000 Subject: Rewrite two abs masks as literals These were technically undefined behaviour - they have been rewritten without the shift so that their type is unsigned int by default. --- pl/math/log10f.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'pl/math/log10f.c') diff --git a/pl/math/log10f.c b/pl/math/log10f.c index 32de42f..5813982 100644 --- a/pl/math/log10f.c +++ b/pl/math/log10f.c @@ -67,7 +67,7 @@ log10f (float x) tmp = ix - OFF; i = (tmp >> (23 - LOGF_TABLE_BITS)) % N; k = (int32_t) tmp >> 23; /* arithmetic shift. */ - iz = ix - (tmp & 0x1ff << 23); + iz = ix - (tmp & 0xff800000); invc = T[i].invc; logc = T[i].logc; z = (double_t) asfloat (iz); -- cgit v1.2.3 From f0f80b8a19b2593491847ed87456694d789f6f80 Mon Sep 17 00:00:00 2001 From: Joe Ramsay Date: Fri, 6 Jan 2023 09:10:57 +0000 Subject: pl/math: Update copyright years All files in pl/math updated to 2023. --- pl/math/log10f.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'pl/math/log10f.c') diff --git a/pl/math/log10f.c b/pl/math/log10f.c index 5813982..5c80008 100644 --- a/pl/math/log10f.c +++ b/pl/math/log10f.c @@ -1,7 +1,7 @@ /* * Single-precision log10 function. * - * Copyright (c) 2022, Arm Limited. + * Copyright (c) 2022-2023, Arm Limited. * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception */ -- cgit v1.2.3