diff options
Diffstat (limited to 'sysdeps/ieee754/flt-32')
| -rw-r--r-- | sysdeps/ieee754/flt-32/e_gammaf_r.c | 331 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/e_log10f.c | 195 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/math_config.h | 2 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/s_exp10m1f.c | 227 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/s_exp2m1f.c | 194 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/s_expm1f.c | 232 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/s_log10p1f.c | 182 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/s_log1pf.c | 271 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/s_log2p1f.c | 248 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/w_log1pf.c | 1 |
10 files changed, 1430 insertions, 453 deletions
diff --git a/sysdeps/ieee754/flt-32/e_gammaf_r.c b/sysdeps/ieee754/flt-32/e_gammaf_r.c index a9730d61c1..6b1f95d50f 100644 --- a/sysdeps/ieee754/flt-32/e_gammaf_r.c +++ b/sysdeps/ieee754/flt-32/e_gammaf_r.c @@ -1,215 +1,176 @@ -/* Implementation of gamma function according to ISO C. - Copyright (C) 1997-2024 Free Software Foundation, Inc. - This file is part of the GNU C Library. +/* Implementation of the gamma function for binary32. - 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. +Copyright (c) 2023-2024 Alexei Sibidanov. - 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. +The original version of this file was copied from the CORE-MATH +project (file src/binary32/tgamma/tgammaf.c, revision a48e352). - You should have received a copy of the GNU Lesser General Public - License along with the GNU C Library; if not, see - <https://www.gnu.org/licenses/>. */ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. + */ + +/* Changes with respect to the original CORE-MATH code: + - removed the dealing with errno + (this is done in the wrapper math/w_tgammaf_compat.c) + - usage of math_narrow_eval to deal with underflow/overflow + - deal with signgamp + */ #include <math.h> -#include <math-narrow-eval.h> -#include <math_private.h> -#include <fenv_private.h> -#include <math-underflow.h> #include <float.h> +#include <stdint.h> +#include <stddef.h> #include <libm-alias-finite.h> +#include <math-narrow-eval.h> +#include "math_config.h" -/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's - approximation to gamma function. */ +float +__ieee754_gammaf_r (float x, int *signgamp) +{ + /* The wrapper in math/w_tgamma_template.c expects *signgamp to be set to a + non-negative value if the returned value is gamma(x), and to a negative + value if it is -gamma(x). + Since the code here directly computes gamma(x), we set it to 1. + */ + if (signgamp != NULL) + *signgamp = 1; -static const float gamma_coeff[] = + /* List of exceptional cases. Each entry contains the 32-bit encoding u of x, + a binary32 approximation f of gamma(x), and a correction term df. */ + static const struct { - 0x1.555556p-4f, - -0xb.60b61p-12f, - 0x3.403404p-12f, + uint32_t u; + float f, df; + } tb[] = { + { 0x27de86a9u, 0x1.268266p+47f, 0x1p22f }, /* x = 0x1.bd0d52p-48 */ + { 0x27e05475u, 0x1.242422p+47f, 0x1p22f }, /* x = 0x1.c0a8eap-48 */ + { 0xb63befb3u, -0x1.5cb6e4p+18f, 0x1p-7f }, /* x = -0x1.77df66p-19 */ + { 0x3c7bb570u, 0x1.021d9p+6f, 0x1p-19f }, /* x = 0x1.f76aep-7 */ + { 0x41e886d1u, 0x1.33136ap+98f, 0x1p73f }, /* x = 0x1.d10da2p+4 */ + { 0xc067d177u, 0x1.f6850cp-3f, 0x1p-28f }, /* x = -0x1.cfa2eep+1 */ + { 0xbd99da31u, -0x1.befe66p+3, -0x1p-22f }, /* x = -0x1.33b462p-4 */ + { 0xbf54c45au, -0x1.a6b4ecp+2, 0x1p-23f }, /* x = -0x1.a988b4p-1 */ + { 0x41ee77feu, 0x1.d3631cp+101, -0x1p-76f }, /* x = 0x1.dceffcp+4 */ + { 0x3f843a64u, 0x1.f6c638p-1, 0x1p-26f }, /* x = 0x1.0874c8p+0 */ }; -#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0])) - -/* Return gamma (X), for positive X less than 42, in the form R * - 2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to - avoid overflow or underflow in intermediate calculations. */ - -static float -gammaf_positive (float x, int *exp2_adj) -{ - int local_signgam; - if (x < 0.5f) - { - *exp2_adj = 0; - return __ieee754_expf (__ieee754_lgammaf_r (x + 1, &local_signgam)) / x; - } - else if (x <= 1.5f) - { - *exp2_adj = 0; - return __ieee754_expf (__ieee754_lgammaf_r (x, &local_signgam)); - } - else if (x < 2.5f) - { - *exp2_adj = 0; - float x_adj = x - 1; - return (__ieee754_expf (__ieee754_lgammaf_r (x_adj, &local_signgam)) - * x_adj); - } - else - { - float eps = 0; - float x_eps = 0; - float x_adj = x; - float prod = 1; - if (x < 4.0f) - { - /* Adjust into the range for applying Stirling's - approximation. */ - float n = ceilf (4.0f - x); - x_adj = math_narrow_eval (x + n); - x_eps = (x - (x_adj - n)); - prod = __gamma_productf (x_adj - n, x_eps, n, &eps); + uint32_t t = asuint (x); + uint32_t ax = t << 1; + if (__glibc_unlikely (ax >= (0xffu << 24))) + { /* x=NaN or +/-Inf */ + if (ax == (0xffu << 24)) + { /* x=+/-Inf */ + if (t >> 31) /* x=-Inf */ + return __math_invalidf (x); + return x; /* x=+Inf */ } - /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)). - Compute gamma (X_ADJ + X_EPS) using Stirling's approximation, - starting by computing pow (X_ADJ, X_ADJ) with a power of 2 - factored out. */ - float exp_adj = -eps; - float x_adj_int = roundf (x_adj); - float x_adj_frac = x_adj - x_adj_int; - int x_adj_log2; - float x_adj_mant = __frexpf (x_adj, &x_adj_log2); - if (x_adj_mant < M_SQRT1_2f) + return x + x; /* x=NaN, where x+x ensures the "Invalid operation" + exception is set if x is sNaN */ + } + double z = x; + if (__glibc_unlikely (ax < 0x6d000000u)) + { /* |x| < 0x1p-18 */ + volatile double d = (0x1.fa658c23b1578p-1 - 0x1.d0a118f324b63p-1 * z) + * z - 0x1.2788cfc6fb619p-1; + double f = 1.0 / z + d; + float r = f; + uint64_t rt = asuint64 (f); + if (((rt + 2) & 0xfffffff) < 4) { - x_adj_log2--; - x_adj_mant *= 2.0f; + for (unsigned i = 0; i < sizeof (tb) / sizeof (tb[0]); i++) + if (t == tb[i].u) + return tb[i].f + tb[i].df; } - *exp2_adj = x_adj_log2 * (int) x_adj_int; - float ret = (__ieee754_powf (x_adj_mant, x_adj) - * __ieee754_exp2f (x_adj_log2 * x_adj_frac) - * __ieee754_expf (-x_adj) - * sqrtf (2 * M_PIf / x_adj) - / prod); - exp_adj += x_eps * __ieee754_logf (x_adj); - float bsum = gamma_coeff[NCOEFF - 1]; - float x_adj2 = x_adj * x_adj; - for (size_t i = 1; i <= NCOEFF - 1; i++) - bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i]; - exp_adj += bsum / x_adj; - return ret + ret * __expm1f (exp_adj); + return r; } -} - -float -__ieee754_gammaf_r (float x, int *signgamp) -{ - int32_t hx; - float ret; - - GET_FLOAT_WORD (hx, x); - - if (__glibc_unlikely ((hx & 0x7fffffff) == 0)) + float fx = floorf (x); + if (__glibc_unlikely (x >= 0x1.18522p+5f)) { - /* Return value for x == 0 is Inf with divide by zero exception. */ - *signgamp = 0; - return 1.0 / x; + /* Overflow case. The original CORE-MATH code returns + 0x1p127f * 0x1p127f, but apparently some compilers replace this + by +Inf. */ + return math_narrow_eval (x * 0x1p127f); } - if (__builtin_expect (hx < 0, 0) - && (uint32_t) hx < 0xff800000 && rintf (x) == x) - { - /* Return value for integer x < 0 is NaN with invalid exception. */ - *signgamp = 0; - return (x - x) / (x - x); + /* compute k only after the overflow check, otherwise the case to integer + might overflow */ + int k = fx; + if (__glibc_unlikely (fx == x)) + { /* x is integer */ + if (x == 0.0f) + return 1.0f / x; + if (x < 0.0f) + return __math_invalidf (0.0f); + double t0 = 1, x0 = 1; + for (int i = 1; i < k; i++, x0 += 1.0) + t0 *= x0; + return t0; } - if (__glibc_unlikely (hx == 0xff800000)) - { - /* x == -Inf. According to ISO this is NaN. */ - *signgamp = 0; - return x - x; + if (__glibc_unlikely (x < -42.0f)) + { /* negative non-integer */ + /* For x < -42, x non-integer, |gamma(x)| < 2^-151. */ + static const float sgn[2] = { 0x1p-127f, -0x1p-127f }; + /* Underflows always happens */ + return math_narrow_eval (0x1p-127f * sgn[k & 1]); } - if (__glibc_unlikely ((hx & 0x7f800000) == 0x7f800000)) + /* The array c[] stores a degree-15 polynomial approximation for + gamma(x). */ + static const double c[] = { - /* Positive infinity (return positive infinity) or NaN (return - NaN). */ - *signgamp = 0; - return x + x; - } + 0x1.c9a76be577123p+0, 0x1.8f2754ddcf90dp+0, 0x1.0d1191949419bp+0, + 0x1.e1f42cf0ae4a1p-2, 0x1.82b358a3ab638p-3, 0x1.e1f2b30cd907bp-5, + 0x1.240f6d4071bd8p-6, 0x1.1522c9f3cd012p-8, 0x1.1fd0051a0525bp-10, + 0x1.9808a8b96c37ep-13, 0x1.b3f78e01152b5p-15, 0x1.49c85a7e1fd04p-18, + 0x1.471ca49184475p-19, -0x1.368f0b7ed9e36p-23, 0x1.882222f9049efp-23, + -0x1.a69ed2042842cp-25 + }; - if (x >= 36.0f) - { - /* Overflow. */ - *signgamp = 0; - ret = math_narrow_eval (FLT_MAX * FLT_MAX); - return ret; - } - else + double m = z - 0x1.7p+1; + double i = roundeven (m); + double step = copysign (1.0, i); + double d = m - i, d2 = d * d, d4 = d2 * d2, d8 = d4 * d4; + double f = (c[0] + d * c[1]) + d2 * (c[2] + d * c[3]) + + d4 * ((c[4] + d * c[5]) + d2 * (c[6] + d * c[7])) + + d8 * ((c[8] + d * c[9]) + d2 * (c[10] + d * c[11]) + + d4 * ((c[12] + d * c[13]) + d2 * (c[14] + d * c[15]))); + int jm = fabs (i); + double w = 1; + if (jm) { - SET_RESTORE_ROUNDF (FE_TONEAREST); - if (x > 0.0f) + z -= 0.5 + step * 0.5; + w = z; + for (int j = jm - 1; j; j--) { - *signgamp = 0; - int exp2_adj; - float tret = gammaf_positive (x, &exp2_adj); - ret = __scalbnf (tret, exp2_adj); + z -= step; + w *= z; } - else if (x >= -FLT_EPSILON / 4.0f) - { - *signgamp = 0; - ret = 1.0f / x; - } - else - { - float tx = truncf (x); - *signgamp = (tx == 2.0f * truncf (tx / 2.0f)) ? -1 : 1; - if (x <= -42.0f) - /* Underflow. */ - ret = FLT_MIN * FLT_MIN; - else - { - float frac = tx - x; - if (frac > 0.5f) - frac = 1.0f - frac; - float sinpix = (frac <= 0.25f - ? __sinf (M_PIf * frac) - : __cosf (M_PIf * (0.5f - frac))); - int exp2_adj; - float tret = M_PIf / (-x * sinpix - * gammaf_positive (-x, &exp2_adj)); - ret = __scalbnf (tret, -exp2_adj); - math_check_force_underflow_nonneg (ret); - } - } - ret = math_narrow_eval (ret); - } - if (isinf (ret) && x != 0) - { - if (*signgamp < 0) - { - ret = math_narrow_eval (-copysignf (FLT_MAX, ret) * FLT_MAX); - ret = -ret; - } - else - ret = math_narrow_eval (copysignf (FLT_MAX, ret) * FLT_MAX); - return ret; } - else if (ret == 0) + if (i <= -0.5) + w = 1 / w; + f *= w; + uint64_t rt = asuint64 (f); + float r = f; + /* Deal with exceptional cases. */ + if (__glibc_unlikely (((rt + 2) & 0xfffffff) < 8)) { - if (*signgamp < 0) - { - ret = math_narrow_eval (-copysignf (FLT_MIN, ret) * FLT_MIN); - ret = -ret; - } - else - ret = math_narrow_eval (copysignf (FLT_MIN, ret) * FLT_MIN); - return ret; + for (unsigned j = 0; j < sizeof (tb) / sizeof (tb[0]); j++) + if (t == tb[j].u) + return tb[j].f + tb[j].df; } - else - return ret; + return r; } libm_alias_finite (__ieee754_gammaf_r, __gammaf_r) diff --git a/sysdeps/ieee754/flt-32/e_log10f.c b/sysdeps/ieee754/flt-32/e_log10f.c index 791895e042..03d9e281f3 100644 --- a/sysdeps/ieee754/flt-32/e_log10f.c +++ b/sysdeps/ieee754/flt-32/e_log10f.c @@ -1,54 +1,161 @@ -/* e_log10f.c -- float version of e_log10.c. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ +/* Correctly-rounded radix-10 logarithm function for binary32 value. + +Copyright (c) 2022-2023 Alexei Sibidanov. + +This file is part of the CORE-MATH project +project (file src/binary32/log10/log10f.c, revision bc385c2). + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ #include <math.h> -#include <math_private.h> -#include <fix-int-fp-convert-zero.h> +#include <stdint.h> #include <libm-alias-finite.h> +#include "math_config.h" -static const float -two25 = 3.3554432000e+07, /* 0x4c000000 */ -ivln10 = 4.3429449201e-01, /* 0x3ede5bd9 */ -log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */ -log10_2lo = 7.9034151668e-07; /* 0x355427db */ +static __attribute__ ((noinline)) float +as_special (float x) +{ + uint32_t ux = asuint (x); + if (ux == 0x7f800000u) + return x; /* +inf */ + uint32_t ax = ux << 1; + if (ax == 0u) + /* -0.0 */ + return __math_divzerof (1); + if (ax > 0xff000000u) + return x + x; /* nan */ + return __math_invalidf (x); +} float -__ieee754_log10f(float x) +__ieee754_log10f (float x) { - float y,z; - int32_t i,k,hx; - - GET_FLOAT_WORD(hx,x); - - k=0; - if (hx < 0x00800000) { /* x < 2**-126 */ - if (__builtin_expect((hx&0x7fffffff)==0, 0)) - return -two25/fabsf (x); /* log(+-0)=-inf */ - if (__builtin_expect(hx<0, 0)) - return (x-x)/(x-x); /* log(-#) = NaN */ - k -= 25; x *= two25; /* subnormal number, scale up x */ - GET_FLOAT_WORD(hx,x); + static const double tr[] = + { + 0x1p+0, 0x1.f81f82p-1, 0x1.f07c1fp-1, 0x1.e9131acp-1, + 0x1.e1e1e1ep-1, 0x1.dae6077p-1, 0x1.d41d41dp-1, 0x1.cd85689p-1, + 0x1.c71c71cp-1, 0x1.c0e0704p-1, 0x1.bacf915p-1, 0x1.b4e81b5p-1, + 0x1.af286bdp-1, 0x1.a98ef6p-1, 0x1.a41a41ap-1, 0x1.9ec8e95p-1, + 0x1.999999ap-1, 0x1.948b0fdp-1, 0x1.8f9c19p-1, 0x1.8acb90fp-1, + 0x1.8618618p-1, 0x1.8181818p-1, 0x1.7d05f41p-1, 0x1.78a4c81p-1, + 0x1.745d174p-1, 0x1.702e05cp-1, 0x1.6c16c17p-1, 0x1.6816817p-1, + 0x1.642c859p-1, 0x1.605816p-1, 0x1.5c9882cp-1, 0x1.58ed231p-1, + 0x1.5555555p-1, 0x1.51d07ebp-1, 0x1.4e5e0a7p-1, 0x1.4afd6ap-1, + 0x1.47ae148p-1, 0x1.446f865p-1, 0x1.4141414p-1, 0x1.3e22cbdp-1, + 0x1.3b13b14p-1, 0x1.3813814p-1, 0x1.3521cfbp-1, 0x1.323e34ap-1, + 0x1.2f684bep-1, 0x1.2c9fb4ep-1, 0x1.29e412ap-1, 0x1.27350b9p-1, + 0x1.2492492p-1, 0x1.21fb781p-1, 0x1.1f7047ep-1, 0x1.1cf06aep-1, + 0x1.1a7b961p-1, 0x1.1811812p-1, 0x1.15b1e5fp-1, 0x1.135c811p-1, + 0x1.1111111p-1, 0x1.0ecf56cp-1, 0x1.0c9715p-1, 0x1.0a6810ap-1, + 0x1.0842108p-1, 0x1.0624dd3p-1, 0x1.041041p-1, 0x1.0204081p-1, + 0.5 + }; + static const double tl[] = + { + -0x1.d45fd6237ebe3p-47, 0x1.b947689311b6ep-8, 0x1.b5e909c96d7d5p-7, + 0x1.45f4f59ed2165p-6, 0x1.af5f92cbd8f1ep-6, 0x1.0ba01a606de8cp-5, + 0x1.3ed119b9a2b7bp-5, 0x1.714834298eec2p-5, 0x1.a30a9d98357fbp-5, + 0x1.d41d512670813p-5, 0x1.02428c0f65519p-4, 0x1.1a23444eecc3ep-4, + 0x1.31b30543f4cb4p-4, 0x1.48f3ed39bfd04p-4, 0x1.5fe8049a0e423p-4, + 0x1.769140a6aa008p-4, 0x1.8cf1836c98cb3p-4, 0x1.a30a9d55541a1p-4, + 0x1.b8de4d1ee823ep-4, 0x1.ce6e4202ca2e6p-4, 0x1.e3bc1accace07p-4, + 0x1.f8c9683b5abd4p-4, 0x1.06cbd68ca9a6ep-3, 0x1.11142f19df73p-3, + 0x1.1b3e71fa7a97fp-3, 0x1.254b4d37a46e3p-3, 0x1.2f3b6912cbf07p-3, + 0x1.390f683115886p-3, 0x1.42c7e7fffc5a8p-3, 0x1.4c65808c78d3cp-3, + 0x1.55e8c50751c55p-3, 0x1.5f52445dec3d8p-3, 0x1.68a288c3f12p-3, + 0x1.71da17bdf0d19p-3, 0x1.7af973608afd9p-3, 0x1.84011952a2579p-3, + 0x1.8cf1837a7ea6p-3, 0x1.95cb2891e43d6p-3, 0x1.9e8e7b0f869ep-3, + 0x1.a73beaa5db18dp-3, 0x1.afd3e394558d3p-3, 0x1.b856cf060d9f1p-3, + 0x1.c0c5134de1ffcp-3, 0x1.c91f1371bc99fp-3, 0x1.d1652ffcd3f53p-3, + 0x1.d997c6f635e75p-3, 0x1.e1b733ab90f3bp-3, 0x1.e9c3ceadac856p-3, + 0x1.f1bdeec43a305p-3, 0x1.f9a5e7a5fa3fep-3, 0x1.00be05ac02f2bp-2, + 0x1.04a054d81a2d4p-2, 0x1.087a0835957fbp-2, 0x1.0c4b457099517p-2, + 0x1.101431aa1fe51p-2, 0x1.13d4f08b98dd8p-2, 0x1.178da53edb892p-2, + 0x1.1b3e71e9f9d58p-2, 0x1.1ee777defdeedp-2, 0x1.2288d7b48e23bp-2, + 0x1.2622b0f52e49fp-2, 0x1.29b522a4c6314p-2, 0x1.2d404b0e30f8p-2, + 0x1.30c4478f3fbe5p-2, 0x1.34413509f7915p-2 + }; + static const union + { + float f; + uint32_t u; + } st[] = + { + { 0x1p+0 }, { 0x1.4p+3 }, { 0x1.9p+6 }, { 0x1.f4p+9 }, + { 0x1.388p+13 }, { 0x1.86ap+16 }, { 0x1.e848p+19 }, { 0x1.312dp+23 }, + { 0x1.7d784p+26 }, { 0x1.dcd65p+29 }, { 0x1.2a05f2p+33 }, { 0 }, + { 0 }, { 0 }, { 0 }, { 0 } + }; + static const double b[] = + { + 0x1.bcb7b15c5a2f8p-2, -0x1.bcbb1dbb88ebap-3, 0x1.2871c39d521c6p-3 + }; + static const double c[] = + { + 0x1.bcb7b1526e50ep-2, -0x1.bcb7b1526e53dp-3, 0x1.287a7636f3fa2p-3, + -0x1.bcb7b146a14b3p-4, 0x1.63c627d5219cbp-4, -0x1.2880736c8762dp-4, + 0x1.fc1ecf913961ap-5 + }; + uint32_t ux = asuint (x); + if (__glibc_unlikely (ux < (1 << 23) || ux >= 0x7f800000u)) + { + if (ux == 0 || ux >= 0x7f800000u) + return as_special (x); + /* subnormal */ + int n = __builtin_clz (ux) - 8; + ux <<= n; + ux -= n << 23; + } + unsigned m = ux & ((1 << 23) - 1), j = (m + (1 << (23 - 7))) >> (23 - 6); + double ix = tr[j], l = tl[j]; + int e = ((int) ux >> 23) - 127; + unsigned je = e + 1; + je = (je * 0x4d104d4) >> 28; + if (__glibc_unlikely (ux == st[je].u)) + return je; + + double tz = asdouble (((int64_t) m | ((int64_t) 1023 << 23)) << (52 - 23)); + double z = tz * ix - 1, z2 = z * z; + double r + = ((e * 0x1.34413509f79ffp-2 + l) + z * b[0]) + z2 * (b[1] + z * b[2]); + float ub = r, lb = r + 0x1.b008p-34; + if (__glibc_unlikely (ub != lb)) + { + double f = z + * ((c[0] + z * c[1]) + + z2 + * ((c[2] + z * c[3]) + + z2 * (c[4] + z * c[5] + z2 * c[6]))); + f -= 0x1.0cee0ed4ca7e9p-54 * e; + f += l - tl[0]; + double el = e * 0x1.34413509f7ap-2; + r = el + f; + ub = r; + tz = r; + if (__glibc_unlikely (!((asuint64 (tz) & ((1 << 28) - 1))))) + { + double dr = (el - r) + f; + r += dr * 32; + ub = r; } - if (__builtin_expect(hx >= 0x7f800000, 0)) return x+x; - k += (hx>>23)-127; - i = ((uint32_t)k&0x80000000)>>31; - hx = (hx&0x007fffff)|((0x7f-i)<<23); - y = (float)(k+i); - if (FIX_INT_FP_CONVERT_ZERO && y == 0.0f) - y = 0.0f; - SET_FLOAT_WORD(x,hx); - z = y*log10_2lo + ivln10*__ieee754_logf(x); - return z+y*log10_2hi; + } + return ub; } libm_alias_finite (__ieee754_log10f, __log10f) diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h index 729f22cd4f..dc07ebd459 100644 --- a/sysdeps/ieee754/flt-32/math_config.h +++ b/sysdeps/ieee754/flt-32/math_config.h @@ -166,9 +166,9 @@ extern const struct exp2f_data uint64_t tab[1 << EXP2F_TABLE_BITS]; double shift_scaled; double poly[EXP2F_POLY_ORDER]; - double shift; double invln2_scaled; double poly_scaled[EXP2F_POLY_ORDER]; + double shift; } __exp2f_data attribute_hidden; #define LOGF_TABLE_BITS 4 diff --git a/sysdeps/ieee754/flt-32/s_exp10m1f.c b/sysdeps/ieee754/flt-32/s_exp10m1f.c new file mode 100644 index 0000000000..ea3173a174 --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_exp10m1f.c @@ -0,0 +1,227 @@ +/* Implementation of the exp10m1 function for binary32. + +Copyright (c) 2022-2024 Alexei Sibidanov. Paul Zimmermann. + +The original version of this file was copied from the CORE-MATH +project (file src/binary32/exp10m1/exp10m1f.c, revision c46b85b). + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. + */ + +#include <math.h> +#include "math_config.h" +#include <libm-alias-float.h> + +float +__exp10m1f (float x) +{ + const double iln10h = 0x1.a934f09p+1 * 16; + const double iln10l = 0x1.e68dc57f2496p-29 * 16; + double z = x; + uint32_t ux = asuint (x); + uint32_t ax = ux & (~0u >> 1); + if (__glibc_unlikely (ux > 0xc0f0d2f1u)) + { /* x < -7.52575 */ + if (ax > (0xffu << 23)) + return x + x; /* nan */ + return (ux == 0xff800000) ? -0x1p+0f : -0x1p+0f + 0x1p-26f; + } + else if (__glibc_unlikely (ax > 0x421a209au)) + { /* x > 38.5318 */ + if (ax >= asuint (INFINITY)) + return x + x; /* +Inf or NaN */ + return __math_oflowf (0); + } + else if (__glibc_unlikely (ax < 0x3d89c604u)) + { /* |x| < 0.1549/log(10) */ + double z2 = z * z, r; + if (__glibc_unlikely (ax < 0x3d1622fbu)) + { /* |x| < 8.44e-2/log(10) */ + if (__glibc_unlikely (ax < 0x3c8b76a3u)) + { /* |x| < 3.92e-2/log(10) */ + if (__glibc_unlikely (ax < 0x3bcced04u)) + { /* |x| < 1.44e-2/log(10) */ + if (__glibc_unlikely (ax < 0x3acf33ebu)) + { /* |x| < 3.64e-3/log(10 */ + if (__glibc_unlikely (ax < 0x395a966bu)) + { /* |x| < 4.8e-4/log(10 */ + if (__glibc_unlikely (ax < 0x36fe4a4bu)) + { /* |x| < 1.745e-5/log(10) */ + if (__glibc_unlikely (ax < 0x32407f39u)) + { /* |x| < 2.58e-8/log(10) */ + if (__glibc_unlikely (ax < 0x245e5bd9u)) + { /* |x| < 4.82164e-17 */ + r = 0x1.26bb1bbb55516p+1; + } + else + { + if (__glibc_unlikely (ux == 0x2c994b7bu)) + return 0x1.60f974p-37f - 0x1p-90f; + r = 0x1.26bb1bbb55516p+1 + + z * 0x1.53524c73cea69p+1; + } + } + else + { + if (__glibc_unlikely (ux == 0xb6fa215bu)) + return -0x1.1ff87ep-16 + 0x1p-68; + r = 0x1.26bb1bbb55516p+1 + + z * (0x1.53524c73ea62fp+1 + + z * 0x1.0470591de2c75p+1); + } + } + else + { + static const double cp[] = + { + 0x1.26bb1bbb55515p+1, 0x1.53524c73cea69p+1, + 0x1.0470595038cc2p+1, 0x1.2bd7609fe1561p+0 + }; + r = (cp[0] + z * cp[1]) + + z2 * (cp[2] + z * cp[3]); + } + } + else + { + static const double cp[] = + { + 0x1.26bb1bbb55516p+1, 0x1.53524c73ce6dbp+1, + 0x1.0470591de3024p+1, 0x1.2bd76b79060e6p+0, + 0x1.1429ffd3a963dp-1 + }; + r = (cp[0] + z * cp[1]) + + z2 * (cp[2] + z * (cp[3] + z * cp[4])); + } + } + else + { + static const double cp[] = + { + 0x1.26bb1bbb55516p+1, 0x1.53524c73cea67p+1, + 0x1.0470591dc2953p+1, 0x1.2bd760a004d64p+0, + 0x1.142a85da6f072p-1, 0x1.a7ed70725b00ep-3 + }; + r = (cp[0] + z * cp[1]) + z2 + * ((cp[2] + z * cp[3]) + + z2 * (cp[4] + z * cp[5])); + } + } + else + { + static const double cp[] = + { + 0x1.26bb1bbb55516p+1, 0x1.53524c73ceadep+1, + 0x1.0470591de2bb4p+1, 0x1.2bd76099a9d33p+0, + 0x1.1429ffd829b0bp-1, 0x1.a7f2a6a0f7dc8p-3, + 0x1.16e4dfbce0f56p-4 + }; + r = (cp[0] + z * cp[1]) + + z2 * ((cp[2] + z * cp[3]) + + z2 * (cp[4] + z * (cp[5] + z * cp[6]))); + } + } + else + { + static const double cp[] = + { + 0x1.26bb1bbb55515p+1, 0x1.53524c73cea6ap+1, + 0x1.0470591de476p+1, 0x1.2bd7609fd4ee2p+0, + 0x1.1429ff70a9b48p-1, 0x1.a7ed71259ba5bp-3, + 0x1.16f3004fb3ac1p-4, 0x1.4116b0388aa9fp-6 + }; + r = ((cp[0] + z * cp[1]) + z2 * (cp[2] + z * cp[3])) + + (z2 * z2) * ((cp[4] + z * cp[5]) + z2 * (cp[6] + + z * cp[7])); + } + } + else + { + static const double cp[] = + { + 0x1.26bb1bbb55515p+1, 0x1.53524c73cea42p+1, 0x1.0470591de2d1dp+1, + 0x1.2bd760a010a53p+0, 0x1.1429ffd16170cp-1, 0x1.a7ed6b2a0d97fp-3, + 0x1.16e4e37fa51e4p-4, 0x1.4147fe4c1676fp-6, 0x1.4897c4b3e329ap-8 + }; + r = ((cp[0] + z * cp[1]) + z2 * (cp[2] + z * cp[3])) + + (z2 * z2) * ((cp[4] + z * cp[5]) + + z2 * (cp[6] + z * (cp[7] + z * cp[8]))); + } + r *= z; + return r; + } + else + { + /* -7.52575 < x < -0.1549/log(10) or 0.1549/log(10) < x < 38.5318 */ + static const double tb[] = + { + 0x1p+0, 0x1.0b5586cf9890fp+0, 0x1.172b83c7d517bp+0, + 0x1.2387a6e756238p+0, 0x1.306fe0a31b715p+0, 0x1.3dea64c123422p+0, + 0x1.4bfdad5362a27p+0, 0x1.5ab07dd485429p+0, 0x1.6a09e667f3bcdp+0, + 0x1.7a11473eb0187p+0, 0x1.8ace5422aa0dap+0, 0x1.9c49182a3f09p+0, + 0x1.ae89f995ad3adp+0, 0x1.c199bdd85529cp+0, 0x1.d5818dcfba487p+0, + 0x1.ea4afa2a490dap+0 + }; + static const double c[] = + { + 0x1.62e42fefa398bp-5, 0x1.ebfbdff84555ap-11, 0x1.c6b08d4ad86d3p-17, + 0x1.3b2ad1b1716a2p-23, 0x1.5d7472718ce9dp-30, 0x1.4a1d7f457ac56p-37 + }; + + if (__glibc_unlikely ((ux << 11) == 0)) + { + uint32_t k = (ux >> 21) - 0x1fc; + if (k <= 0xb) + { + if (k == 0) + return 10.0f - 1.0f; + if (k == 4) + return 100.0f - 1.0f; + if (k == 6) + return 1000.0f - 1.0f; + if (k == 8) + return 10000.0f - 1.0f; + if (k == 9) + return 100000.0f - 1.0f; + if (k == 10) + return 1000000.0f - 1.0f; + if (k == 11) + return 10000000.0f - 1.0f; + } + } + double a = iln10h * z; + double ia = floor (a); + double h = (a - ia) + iln10l * z; + int64_t i = ia; + int64_t j = i & 0xf; + int64_t e = i - j; + e >>= 4; + double s = tb[j]; + s *= asdouble ((e + 0x3ffull) << 52); + double h2 = h * h; + double c0 = c[0] + h * c[1]; + double c2 = c[2] + h * c[3]; + double c4 = c[4] + h * c[5]; + c0 += h2 * (c2 + h2 * c4); + double w = s * h; + return (s - 1.0) + w * c0; + } +} +#ifndef __exp10m1f +libm_alias_float (__exp10m1, exp10m1) +#endif diff --git a/sysdeps/ieee754/flt-32/s_exp2m1f.c b/sysdeps/ieee754/flt-32/s_exp2m1f.c new file mode 100644 index 0000000000..325ffb11b0 --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_exp2m1f.c @@ -0,0 +1,194 @@ +/* Correctly-rounded base-2 exponent function biased by 1 for binary32 value. + +Copyright (c) 2022-2024 Alexei Sibidanov. + +The original version of this file was copied from the CORE-MATH +project (file src/binary32/exp2m1/exp2m1f.c, revision baf5f6b). + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ + +#include <fenv.h> +#include <math.h> +#include "math_config.h" +#include <libm-alias-float.h> +#include <math-narrow-eval.h> +#include <float.h> + +float +__exp2m1f (float x) +{ + double z = x; + uint32_t ux = asuint (x); + uint32_t ax = ux & (~0u >> 1); + if (__glibc_unlikely (ux >= 0xc1c80000u)) + { /* x <= -25 */ + if (ax > (0xffu << 23)) + return x + x; /* nan */ + return (ux == 0xff800000) ? -0x1p+0f : -0x1p+0f + 0x1p-26f; + } + else if (__glibc_unlikely (ax >= 0x43000000u)) + { /* x >= 128 */ + if (ax >= asuint (INFINITY)) + return x + x; /* +Inf or NaN */ + /* exp2m1 (MAX_EXP) should not overflow when rounding towards zero + or towards -Inf. We round FLT_MAX + 2^103 which is in the middle + between FLT_MAX and 2^128 (the next number with unbounded range). */ + float ret = math_narrow_eval (FLT_MAX + 0x1p103f); + if (x == FLT_MAX_EXP && ret == FLT_MAX) + return ret; + return __math_oflowf (0); + } + else if (__glibc_unlikely (ax < 0x3df95f1fu)) + { /* |x| < 8.44e-2/log(2) */ + double z2 = z * z, r; + if (__glibc_unlikely (ax < 0x3d67a4ccu)) + { /* |x| < 3.92e-2/log(2) */ + if (__glibc_unlikely (ax < 0x3caa2feeu)) + { /* |x| < 1.44e-2/log(2) */ + if (__glibc_unlikely (ax < 0x3bac1405u)) + { /* |x| < 3.64e-3/log(2) */ + if (__glibc_unlikely (ax < 0x3a358876u)) + { /* |x| < 4.8e-4/log(2) */ + if (__glibc_unlikely (ax < 0x37d32ef6u)) + { /* |x| < 1.745e-5/log(2) */ + if (__glibc_unlikely (ax < 0x331fdd82u)) + { /* |x| < 2.58e-8/log(2) */ + if (__glibc_unlikely (ax < 0x2538aa3bu)) + /* |x| < 1.60171e-16 */ + r = 0x1.62e42fefa39efp-1; + else + r = 0x1.62e42fefa39fp-1 + + z * 0x1.ebfbdff82c58fp-3; + } + else + { + if (__glibc_unlikely (ux == 0xb3d85005u)) + return -0x1.2bdf76p-24 - 0x1.8p-77; + if (__glibc_unlikely (ux == 0x3338428du)) + return 0x1.fee08ap-26 + 0x1p-80; + static const double c[] = + { + 0x1.62e42fefa39efp-1, 0x1.ebfbdff8548fdp-3, + 0x1.c6b08d704a06dp-5 + }; + r = c[0] + z * (c[1] + z * c[2]); + } + } + else + { + if (__glibc_unlikely (ux == 0x388bca4fu)) + return 0x1.839702p-15 - 0x1.8p-68; + static const double c[] = + { + 0x1.62e42fefa39efp-1, 0x1.ebfbdff82c58fp-3, + 0x1.c6b08dc82b347p-5, 0x1.3b2ab6fbad172p-7 + }; + r = (c[0] + z * c[1]) + z2 * (c[2] + z * c[3]); + } + } + else + { + static const double c[] = + { + 0x1.62e42fefa39efp-1, 0x1.ebfbdff82c068p-3, + 0x1.c6b08d704a6dcp-5, 0x1.3b2ac262c3eedp-7, + 0x1.5d87fe7af779ap-10 + }; + r = (c[0] + z * c[1]) + + z2 * (c[2] + z * (c[3] + z * c[4])); + } + } + else + { + static const double c[] = + { + 0x1.62e42fefa39fp-1, 0x1.ebfbdff82c58dp-3, + 0x1.c6b08d7011d13p-5, 0x1.3b2ab6fbd267dp-7, + 0x1.5d88a81cea49ep-10, 0x1.430912ea9b963p-13 + }; + r = (c[0] + z * c[1]) + + z2 * ((c[2] + z * c[3]) + z2 * (c[4] + z * c[5])); + } + } + else + { + static const double c[] = + { + 0x1.62e42fefa39efp-1, 0x1.ebfbdff82c639p-3, + 0x1.c6b08d7049f1cp-5, 0x1.3b2ab6f5243bdp-7, + 0x1.5d87fe80a9e6cp-10, 0x1.430d0b9257fa8p-13, + 0x1.ffcbfc4cf0952p-17 + }; + r = (c[0] + z * c[1]) + + z2 * ((c[2] + z * c[3]) + + z2 * (c[4] + z * (c[5] + z * c[6]))); + } + } + else + { + static const double c[] = + { + 0x1.62e42fefa39efp-1, 0x1.ebfbdff82c591p-3, + 0x1.c6b08d704cf6bp-5, 0x1.3b2ab6fba00cep-7, + 0x1.5d87fdfdaadb4p-10, 0x1.4309137333066p-13, + 0x1.ffe5e90daf7ddp-17, 0x1.62c0220eed731p-20 + }; + r = ((c[0] + z * c[1]) + z2 * (c[2] + z * c[3])) + + (z2 * z2) * ((c[4] + z * c[5]) + z2 * (c[6] + z * c[7])); + } + r *= z; + return r; + } + else + { + static const double c[] = + { + 0x1.62e42fefa398bp-5, 0x1.ebfbdff84555ap-11, + 0x1.c6b08d4ad86d3p-17, 0x1.3b2ad1b1716a2p-23, + 0x1.5d7472718ce9dp-30, 0x1.4a1d7f457ac56p-37 + }; + static const double tb[] = + { + 0x1p+0, 0x1.0b5586cf9890fp+0, 0x1.172b83c7d517bp+0, + 0x1.2387a6e756238p+0, 0x1.306fe0a31b715p+0, 0x1.3dea64c123422p+0, + 0x1.4bfdad5362a27p+0, 0x1.5ab07dd485429p+0, 0x1.6a09e667f3bcdp+0, + 0x1.7a11473eb0187p+0, 0x1.8ace5422aa0dap+0, 0x1.9c49182a3f09p+0, + 0x1.ae89f995ad3adp+0, 0x1.c199bdd85529cp+0, 0x1.d5818dcfba487p+0, + 0x1.ea4afa2a490dap+0 + }; + double a = 16.0 * z; + double ia = floor (a); + double h = a - ia; + double h2 = h * h; + int64_t i = ia, j = i & 0xf, e = i - j; + e >>= 4; + double s = tb[j]; + s *= asdouble ((e + 0x3ffull) << 52); + double c0 = c[0] + h * c[1]; + double c2 = c[2] + h * c[3]; + double c4 = c[4] + h * c[5]; + c0 += h2 * (c2 + h2 * c4); + double w = s * h; + return (s - 1.0) + w * c0; + } +} +#ifndef __exp2m1f +libm_alias_float (__exp2m1, exp2m1) +#endif diff --git a/sysdeps/ieee754/flt-32/s_expm1f.c b/sysdeps/ieee754/flt-32/s_expm1f.c index 35f7b5214a..edd7c9acf8 100644 --- a/sysdeps/ieee754/flt-32/s_expm1f.c +++ b/sysdeps/ieee754/flt-32/s_expm1f.c @@ -1,132 +1,124 @@ -/* s_expm1f.c -- float version of s_expm1.c. - */ +/* Correctly-rounded natural exponent function biased by 1 for binary32 value. -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ +Copyright (c) 2022-2024 Alexei Sibidanov. + +This file is part of the CORE-MATH project +project (file src/binary32/expm1/expm1f.c, revision bc385c2). + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ -#include <errno.h> -#include <float.h> #include <math.h> -#include <math-barriers.h> -#include <math_private.h> #include <math-underflow.h> #include <libm-alias-float.h> - -static const float huge = 1.0e+30; -static const float tiny = 1.0e-30; - -static const float -one = 1.0, -o_threshold = 8.8721679688e+01,/* 0x42b17180 */ -ln2_hi = 6.9313812256e-01,/* 0x3f317180 */ -ln2_lo = 9.0580006145e-06,/* 0x3717f7d1 */ -invln2 = 1.4426950216e+00,/* 0x3fb8aa3b */ - /* scaled coefficients related to expm1 */ -Q1 = -3.3333335072e-02, /* 0xbd088889 */ -Q2 = 1.5873016091e-03, /* 0x3ad00d01 */ -Q3 = -7.9365076090e-05, /* 0xb8a670cd */ -Q4 = 4.0082177293e-06, /* 0x36867e54 */ -Q5 = -2.0109921195e-07; /* 0xb457edbb */ +#include "math_config.h" float -__expm1f(float x) +__expm1f (float x) { - float y,hi,lo,c,t,e,hxs,hfx,r1; - int32_t k,xsb; - uint32_t hx; - - GET_FLOAT_WORD(hx,x); - xsb = hx&0x80000000; /* sign bit of x */ - if(xsb==0) y=x; else y= -x; /* y = |x| */ - hx &= 0x7fffffff; /* high word of |x| */ - - /* filter out huge and non-finite argument */ - if(hx >= 0x4195b844) { /* if |x|>=27*ln2 */ - if(hx >= 0x42b17218) { /* if |x|>=88.721... */ - if(hx>0x7f800000) - return x+x; /* NaN */ - if(hx==0x7f800000) - return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */ - if(x > o_threshold) { - __set_errno (ERANGE); - return huge*huge; /* overflow */ - } - } - if(xsb!=0) { /* x < -27*ln2, return -1.0 with inexact */ - math_force_eval(x+tiny);/* raise inexact */ - return tiny-one; /* return -1 */ - } - } - - /* argument reduction */ - if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ - if(hx < 0x3F851592) { /* and |x| < 1.5 ln2 */ - if(xsb==0) - {hi = x - ln2_hi; lo = ln2_lo; k = 1;} - else - {hi = x + ln2_hi; lo = -ln2_lo; k = -1;} - } else { - k = invln2*x+((xsb==0)?(float)0.5:(float)-0.5); - t = k; - hi = x - t*ln2_hi; /* t*ln2_hi is exact here */ - lo = t*ln2_lo; - } - x = hi - lo; - c = (hi-x)-lo; - } - else if(hx < 0x33000000) { /* when |x|<2**-25, return x */ - math_check_force_underflow (x); - t = huge+x; /* return x with inexact flags when x!=0 */ - return x - (t-(huge+x)); + static const double c[] = + { + 1, 0x1.62e42fef4c4e7p-6, 0x1.ebfd1b232f475p-13, 0x1.c6b19384ecd93p-20 + }; + static const double ch[] = + { + 0x1.62e42fefa39efp-6, 0x1.ebfbdff82c58fp-13, 0x1.c6b08d702e0edp-20, + 0x1.3b2ab6fb92e5ep-27, 0x1.5d886e6d54203p-35, 0x1.430976b8ce6efp-43 + }; + static const double td[] = + { + 0x1p+0, 0x1.059b0d3158574p+0, 0x1.0b5586cf9890fp+0, + 0x1.11301d0125b51p+0, 0x1.172b83c7d517bp+0, 0x1.1d4873168b9aap+0, + 0x1.2387a6e756238p+0, 0x1.29e9df51fdee1p+0, 0x1.306fe0a31b715p+0, + 0x1.371a7373aa9cbp+0, 0x1.3dea64c123422p+0, 0x1.44e086061892dp+0, + 0x1.4bfdad5362a27p+0, 0x1.5342b569d4f82p+0, 0x1.5ab07dd485429p+0, + 0x1.6247eb03a5585p+0, 0x1.6a09e667f3bcdp+0, 0x1.71f75e8ec5f74p+0, + 0x1.7a11473eb0187p+0, 0x1.82589994cce13p+0, 0x1.8ace5422aa0dbp+0, + 0x1.93737b0cdc5e5p+0, 0x1.9c49182a3f09p+0, 0x1.a5503b23e255dp+0, + 0x1.ae89f995ad3adp+0, 0x1.b7f76f2fb5e47p+0, 0x1.c199bdd85529cp+0, + 0x1.cb720dcef9069p+0, 0x1.d5818dcfba487p+0, 0x1.dfc97337b9b5fp+0, + 0x1.ea4afa2a490dap+0, 0x1.f50765b6e454p+0 + }; + const double iln2 = 0x1.71547652b82fep+5; + const double big = 0x1.8p52; + double z = x; + uint32_t ux = asuint (x); + uint32_t ax = ux << 1; + if (__glibc_likely (ax < 0x7c400000u)) + { /* |x| < 0.15625 */ + if (__glibc_unlikely (ax < 0x676a09e8u)) + { /* |x| < 0x1.6a09e8p-24 */ + if (__glibc_unlikely (ax == 0x0u)) + return x; /* x = +-0 */ + return fmaf (fabsf (x), 0x1p-25f, x); } - else k = 0; - - /* x is now in primary range */ - hfx = (float)0.5*x; - hxs = x*hfx; - r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5)))); - t = (float)3.0-r1*hfx; - e = hxs*((r1-t)/((float)6.0 - x*t)); - if(k==0) return x - (x*e-hxs); /* c is 0 */ - else { - e = (x*(e-c)-c); - e -= hxs; - if(k== -1) return (float)0.5*(x-e)-(float)0.5; - if(k==1) { - if(x < (float)-0.25) return -(float)2.0*(e-(x+(float)0.5)); - else return one+(float)2.0*(x-e); - } - if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */ - int32_t i; - y = one-(e-x); - GET_FLOAT_WORD(i,y); - SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ - return y-one; - } - t = one; - if(k<23) { - int32_t i; - SET_FLOAT_WORD(t,0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */ - y = t-(e-x); - GET_FLOAT_WORD(i,y); - SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ - } else { - int32_t i; - SET_FLOAT_WORD(t,((0x7f-k)<<23)); /* 2^-k */ - y = x-(e+t); - y += one; - GET_FLOAT_WORD(i,y); - SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ - } + static const double b[] = + { + 0x1.fffffffffffc2p-2, 0x1.55555555555fep-3, 0x1.555555559767fp-5, + 0x1.1111111098dc1p-7, 0x1.6c16bca988aa9p-10, 0x1.a01a07658483fp-13, + 0x1.a05b04d2c3503p-16, 0x1.71de3a960b5e3p-19 + }; + double z2 = z * z, z4 = z2 * z2; + double r = z + z2 + * ((b[0] + z * b[1]) + z2 * (b[2] + z * b[3]) + + z4 * ((b[4] + z * b[5]) + z2 * (b[6] + z * b[7]))); + return r; + } + if (__glibc_unlikely (ax >= 0x8562e430u)) + { /* |x| > 88.72 */ + if (ax > (0xffu << 24)) + return x + x; /* nan */ + if (__glibc_unlikely (ux >> 31)) + { /* x < 0 */ + if (ax == (0xffu << 24)) + return -1.0f; + return -1.0f + 0x1p-26f; } - return y; + if (ax == (0xffu << 24)) + return INFINITY; + return __math_oflowf (0); + } + double a = iln2 * z; + double ia = roundeven (a); + double h = a - ia; + double h2 = h * h; + uint64_t u = asuint64 (ia + big); + double c2 = c[2] + h * c[3], c0 = c[0] + h * c[1]; + const uint64_t *tdl = (uint64_t *) ((void *) td); + double sv = asdouble (tdl[u & 0x1f] + ((u >> 5) << 52)); + double r = (c0 + h2 * c2) * sv - 1.0; + float ub = r, lb = r - sv * 0x1.3b3p-33; + if (__glibc_unlikely (ub != lb)) + { + if (__glibc_unlikely (ux > 0xc18aa123u)) /* x < -17.32 */ + return -1.0f + 0x1p-26f; + const double iln2h = 0x1.7154765p+5; + const double iln2l = 0x1.5c17f0bbbe88p-26; + double s = sv; + h = (iln2h * z - ia) + iln2l * z; + h2 = h * h; + double w = s * h; + r = (s - 1) + w + * ((ch[0] + h * ch[1]) + + h2 * ((ch[2] + h * ch[3]) + h2 * (ch[4] + h * ch[5]))); + ub = r; + } + return ub; } libm_alias_float (__expm1, expm1) diff --git a/sysdeps/ieee754/flt-32/s_log10p1f.c b/sysdeps/ieee754/flt-32/s_log10p1f.c new file mode 100644 index 0000000000..64deb1eeda --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_log10p1f.c @@ -0,0 +1,182 @@ +/* Correctly-rounded biased argument base-10 logarithm function for binary32 value. + +Copyright (c) 2022-2023 Alexei Sibidanov. + +This file is part of the CORE-MATH project +project (file src/binary32/log10p1/log10p1f.c revision bc385c2). + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ + +#include <stdint.h> +#include <errno.h> +#include <fenv.h> +#include <libm-alias-float.h> +#include "math_config.h" + +static __attribute__ ((noinline)) float +as_special (float x) +{ + uint32_t ux = asuint (x); + if (ux == 0x7f800000u) + return x; /* +inf */ + uint32_t ax = ux << 1; + if (ax == 0x17fu << 24) + /* x+1 = 0.0 */ + return __math_divzerof (1); + if (ax > 0xff000000u) + return x + x; /* nan */ + return __math_invalidf (x); +} + +float +__log10p1f (float x) +{ + static const double tr[] = + { + 0x1p+0, 0x1.f81f82p-1, 0x1.f07c1fp-1, 0x1.e9131acp-1, + 0x1.e1e1e1ep-1, 0x1.dae6077p-1, 0x1.d41d41dp-1, 0x1.cd85689p-1, + 0x1.c71c71cp-1, 0x1.c0e0704p-1, 0x1.bacf915p-1, 0x1.b4e81b5p-1, + 0x1.af286bdp-1, 0x1.a98ef6p-1, 0x1.a41a41ap-1, 0x1.9ec8e95p-1, + 0x1.999999ap-1, 0x1.948b0fdp-1, 0x1.8f9c19p-1, 0x1.8acb90fp-1, + 0x1.8618618p-1, 0x1.8181818p-1, 0x1.7d05f41p-1, 0x1.78a4c81p-1, + 0x1.745d174p-1, 0x1.702e05cp-1, 0x1.6c16c17p-1, 0x1.6816817p-1, + 0x1.642c859p-1, 0x1.605816p-1, 0x1.5c9882cp-1, 0x1.58ed231p-1, + 0x1.5555555p-1, 0x1.51d07ebp-1, 0x1.4e5e0a7p-1, 0x1.4afd6ap-1, + 0x1.47ae148p-1, 0x1.446f865p-1, 0x1.4141414p-1, 0x1.3e22cbdp-1, + 0x1.3b13b14p-1, 0x1.3813814p-1, 0x1.3521cfbp-1, 0x1.323e34ap-1, + 0x1.2f684bep-1, 0x1.2c9fb4ep-1, 0x1.29e412ap-1, 0x1.27350b9p-1, + 0x1.2492492p-1, 0x1.21fb781p-1, 0x1.1f7047ep-1, 0x1.1cf06aep-1, + 0x1.1a7b961p-1, 0x1.1811812p-1, 0x1.15b1e5fp-1, 0x1.135c811p-1, + 0x1.1111111p-1, 0x1.0ecf56cp-1, 0x1.0c9715p-1, 0x1.0a6810ap-1, + 0x1.0842108p-1, 0x1.0624dd3p-1, 0x1.041041p-1, 0x1.0204081p-1, + 0.5 + }; + static const double tl[] = + { + 0x1.562ec497ef351p-43, 0x1.b9476892ea99cp-8, 0x1.b5e909c959eecp-7, + 0x1.45f4f59ec84fp-6, 0x1.af5f92cbcf2aap-6, 0x1.0ba01a6069052p-5, + 0x1.3ed119b99dd41p-5, 0x1.714834298a088p-5, 0x1.a30a9d98309c1p-5, + 0x1.d41d51266b9d9p-5, 0x1.02428c0f62dfcp-4, 0x1.1a23444eea521p-4, + 0x1.31b30543f2597p-4, 0x1.48f3ed39bd5e7p-4, 0x1.5fe8049a0bd06p-4, + 0x1.769140a6a78eap-4, 0x1.8cf1836c96595p-4, 0x1.a30a9d5551a84p-4, + 0x1.b8de4d1ee5b21p-4, 0x1.ce6e4202c7bc9p-4, 0x1.e3bc1accaa6eap-4, + 0x1.f8c9683b584b7p-4, 0x1.06cbd68ca86ep-3, 0x1.11142f19de3a2p-3, + 0x1.1b3e71fa795fp-3, 0x1.254b4d37a3354p-3, 0x1.2f3b6912cab79p-3, + 0x1.390f6831144f7p-3, 0x1.42c7e7fffb21ap-3, 0x1.4c65808c779aep-3, + 0x1.55e8c507508c7p-3, 0x1.5f52445deb049p-3, 0x1.68a288c3efe72p-3, + 0x1.71da17bdef98bp-3, 0x1.7af9736089c4bp-3, 0x1.84011952a11ebp-3, + 0x1.8cf1837a7d6d1p-3, 0x1.95cb2891e3048p-3, 0x1.9e8e7b0f85651p-3, + 0x1.a73beaa5d9dfep-3, 0x1.afd3e39454544p-3, 0x1.b856cf060c662p-3, + 0x1.c0c5134de0c6dp-3, 0x1.c91f1371bb611p-3, 0x1.d1652ffcd2bc5p-3, + 0x1.d997c6f634ae6p-3, 0x1.e1b733ab8fbadp-3, 0x1.e9c3ceadab4c8p-3, + 0x1.f1bdeec438f77p-3, 0x1.f9a5e7a5f906fp-3, 0x1.00be05ac02564p-2, + 0x1.04a054d81990cp-2, 0x1.087a083594e33p-2, 0x1.0c4b457098b4fp-2, + 0x1.101431aa1f48ap-2, 0x1.13d4f08b98411p-2, 0x1.178da53edaecbp-2, + 0x1.1b3e71e9f9391p-2, 0x1.1ee777defd526p-2, 0x1.2288d7b48d874p-2, + 0x1.2622b0f52dad8p-2, 0x1.29b522a4c594cp-2, 0x1.2d404b0e305b9p-2, + 0x1.30c4478f3f21dp-2, 0x1.34413509f6f4dp-2 + }; + static const union + { + float f; + uint32_t u; + } st[] = + { + { 0x0p+0 }, { 0x1.2p+3 }, { 0x1.8cp+6 }, + { 0x1.f38p+9 }, { 0x1.3878p+13 }, { 0x1.869fp+16 }, + { 0x1.e847ep+19 }, { 0x1.312cfep+23 } + }; + double z = x; + uint32_t ux = asuint (x); + if (__glibc_unlikely (ux >= 0x17fu << 23)) /* x <= -1 */ + return as_special (x); + uint32_t ax = ux & (~0u >> 1); + if (__glibc_unlikely (ax == 0)) + return copysign (0, x); + if (__glibc_unlikely (ax >= (0xff << 23))) /* +inf, nan */ + return as_special (x); + int ie = ux; + ie >>= 23; + unsigned int je = ie - 126; + je = (je * 0x9a209a8) >> 29; + if (__glibc_unlikely (ux == st[je].u)) + return je; + + uint64_t tz = asuint64 (z + 1.0); + uint64_t m = tz & (~(uint64_t) 0 >> 12); + int32_t e = (tz >> 52) - 1023, j = ((m + ((int64_t) 1 << 45)) >> 46); + tz = m | ((uint64_t) 0x3ff << 52); + double ix = tr[j], l = tl[j]; + double off = e * 0x1.34413509f79ffp-2 + l; + double v = asdouble (tz) * ix - 1; + + static const double h[] = + { + 0x1.bcb7b150bf6d8p-2, -0x1.bcb7b1738c07ep-3, + 0x1.287de19e795c5p-3, -0x1.bca44edc44bc4p-4 + }; + double v2 = v * v; + double f = (h[0] + v * h[1]) + v2 * (h[2] + v * h[3]); + double r = off + v * f; + float ub = r; + float lb = r + 0x1.5cp-42; + if (__glibc_unlikely (ub != lb)) + { + if (__glibc_unlikely (ax < 0x3d32743eu)) + { /* 0x1.64e87cp-5f */ + if (__glibc_unlikely (ux == 0xa6aba8afu)) + return -0x1.2a33bcp-51f + 0x1p-76f; + if (__glibc_unlikely (ux == 0xaf39b9a7u)) + return -0x1.42a342p-34f + 0x1p-59f; + if (__glibc_unlikely (ux == 0x399a7c00u)) + return 0x1.0c53cap-13f + 0x1p-38f; + z /= 2.0 + z; + double z2 = z * z, z4 = z2 * z2; + static const double c[] = + { + 0x1.bcb7b1526e50fp-1, 0x1.287a76370129dp-2, + 0x1.63c62378fa3dbp-3, 0x1.fca4139a42374p-4 + }; + float ret = z * ((c[0] + z2 * c[1]) + z4 * (c[2] + z2 * c[3])); + if (x != 0.0f && ret == 0.0) + __set_errno (ERANGE); + return ret; + } + if (__glibc_unlikely (ux == 0x7956ba5eu)) + return 0x1.16bebap+5f + 0x1p-20f; + if (__glibc_unlikely (ux == 0xbd86ffb9u)) + return -0x1.e53536p-6f + 0x1p-31f; + static const double c[] = + { + 0x1.bcb7b1526e50ep-2, -0x1.bcb7b1526e53dp-3, 0x1.287a7636f3fa2p-3, + -0x1.bcb7b146a14b3p-4, 0x1.63c627d5219cbp-4, -0x1.2880736c8762dp-4, + 0x1.fc1ecf913961ap-5 + }; + f = v + * ((c[0] + v * c[1]) + + v2 * ((c[2] + v * c[3]) + v2 * (c[4] + v * c[5] + v2 * c[6]))); + f += l - tl[0]; + double el = e * 0x1.34413509f79ffp-2; + r = el + f; + ub = r; + } + return ub; +} +libm_alias_float (__log10p1, log10p1) diff --git a/sysdeps/ieee754/flt-32/s_log1pf.c b/sysdeps/ieee754/flt-32/s_log1pf.c index 555f0f82c8..d1686e78aa 100644 --- a/sysdeps/ieee754/flt-32/s_log1pf.c +++ b/sysdeps/ieee754/flt-32/s_log1pf.c @@ -1,116 +1,181 @@ -/* s_log1pf.c -- float version of s_log1p.c. - */ +/* Correctly-rounded biased argument natural logarithm function for binary32 + value. -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ +Copyright (c) 2023, 2024 Alexei Sibidanov. -#include <float.h> -#include <math.h> -#include <math-barriers.h> -#include <math_private.h> -#include <math-underflow.h> -#include <libc-diag.h> +This file is part of the CORE-MATH project +project (file src/binary32/log1p/log1pf.c revision bc385c2). -static const float -ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ -ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ -two25 = 3.355443200e+07, /* 0x4c000000 */ -Lp1 = 6.6666668653e-01, /* 3F2AAAAB */ -Lp2 = 4.0000000596e-01, /* 3ECCCCCD */ -Lp3 = 2.8571429849e-01, /* 3E924925 */ -Lp4 = 2.2222198546e-01, /* 3E638E29 */ -Lp5 = 1.8183572590e-01, /* 3E3A3325 */ -Lp6 = 1.5313838422e-01, /* 3E1CD04F */ -Lp7 = 1.4798198640e-01; /* 3E178897 */ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: -static const float zero = 0.0; +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. -float -__log1pf(float x) +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ + +#include <math.h> +#include <stdint.h> +#include <errno.h> +#include <libm-alias-float.h> +#include "math_config.h" + +static __attribute__ ((noinline)) float +as_special (float x) { - float hfsq,f,c,s,z,R,u; - int32_t k,hx,hu,ax; + uint32_t t = asuint (x); + if (t == 0xbf800000u) + return __math_divzerof (1); + if (t == 0x7f800000u) + return x; /* +inf */ + uint32_t ax = t << 1; + if (ax > 0xff000000u) + return x + x; /* nan */ + return __math_invalidf (0.0f); +} - GET_FLOAT_WORD(hx,x); - ax = hx&0x7fffffff; +float +__log1pf (float x) +{ + static const double x0[] = + { + 0x1.f81f82p-1, 0x1.e9131acp-1, 0x1.dae6077p-1, 0x1.cd85689p-1, + 0x1.c0e0704p-1, 0x1.b4e81b5p-1, 0x1.a98ef6p-1, 0x1.9ec8e95p-1, + 0x1.948b0fdp-1, 0x1.8acb90fp-1, 0x1.8181818p-1, 0x1.78a4c81p-1, + 0x1.702e05cp-1, 0x1.6816817p-1, 0x1.605816p-1, 0x1.58ed231p-1, + 0x1.51d07ebp-1, 0x1.4afd6ap-1, 0x1.446f865p-1, 0x1.3e22cbdp-1, + 0x1.3813814p-1, 0x1.323e34ap-1, 0x1.2c9fb4ep-1, 0x1.27350b9p-1, + 0x1.21fb781p-1, 0x1.1cf06aep-1, 0x1.1811812p-1, 0x1.135c811p-1, + 0x1.0ecf56cp-1, 0x1.0a6810ap-1, 0x1.0624dd3p-1, 0x1.0204081p-1 + }; + static const double lixb[] = + { + 0x1.fc0a8909b4218p-7, 0x1.77458f51aac89p-5, 0x1.341d793afb997p-4, + 0x1.a926d3a5ebd2ap-4, 0x1.0d77e7a8a823dp-3, 0x1.44d2b6c557102p-3, + 0x1.7ab89040accecp-3, 0x1.af3c94ecab3d6p-3, 0x1.e27076d54e6c9p-3, + 0x1.0a324e3888ad5p-2, 0x1.22941fc0c7357p-2, 0x1.3a64c56ae3fdbp-2, + 0x1.51aad874af21fp-2, 0x1.686c81d300eap-2, 0x1.7eaf83c7fa9b5p-2, + 0x1.947941aa610ecp-2, 0x1.a9cec9a3f023bp-2, 0x1.beb4d9ea4156ep-2, + 0x1.d32fe7f35e5c7p-2, 0x1.e7442617b817ap-2, 0x1.faf588dd5ed1p-2, + 0x1.0723e5c635c39p-1, 0x1.109f39d53c99p-1, 0x1.19ee6b38a4668p-1, + 0x1.23130d7f93c3bp-1, 0x1.2c0e9ec9b0b85p-1, 0x1.34e289cb35eccp-1, + 0x1.3d9026ad3d3f3p-1, 0x1.4618bc1eadbbbp-1, 0x1.4e7d8127dd8a9p-1, + 0x1.56bf9d5967092p-1, 0x1.5ee02a926936ep-1 + }; + static const double lix[] = + { + 0x1.fc0a890fc03e4p-7, 0x1.77458f532dcfcp-5, 0x1.341d793bbd1d1p-4, + 0x1.a926d3a6ad563p-4, 0x1.0d77e7a908e59p-3, 0x1.44d2b6c5b7d1ep-3, + 0x1.7ab890410d909p-3, 0x1.af3c94ed0bff3p-3, 0x1.e27076d5af2e6p-3, + 0x1.0a324e38b90e3p-2, 0x1.22941fc0f7966p-2, 0x1.3a64c56b145eap-2, + 0x1.51aad874df82dp-2, 0x1.686c81d3314afp-2, 0x1.7eaf83c82afc3p-2, + 0x1.947941aa916fbp-2, 0x1.a9cec9a42084ap-2, 0x1.beb4d9ea71b7cp-2, + 0x1.d32fe7f38ebd5p-2, 0x1.e7442617e8788p-2, 0x1.faf588dd8f31fp-2, + 0x1.0723e5c64df4p-1, 0x1.109f39d554c97p-1, 0x1.19ee6b38bc96fp-1, + 0x1.23130d7fabf43p-1, 0x1.2c0e9ec9c8e8cp-1, 0x1.34e289cb4e1d3p-1, + 0x1.3d9026ad556fbp-1, 0x1.4618bc1ec5ec2p-1, 0x1.4e7d8127f5bb1p-1, + 0x1.56bf9d597f399p-1, 0x1.5ee02a9281675p-1 + }; + static const double b[] = + { + 0x1p+0, + -0x1p-1, + 0x1.5555555556f6bp-2, + -0x1.00000000029b9p-2, + 0x1.9999988d176e4p-3, + -0x1.55555418889a7p-3, + 0x1.24adeca50e2bcp-3, + -0x1.001ba33bf57cfp-3 + }; - k = 1; - if (hx < 0x3ed413d7) { /* x < 0.41422 */ - if(ax>=0x3f800000) { /* x <= -1.0 */ - if(x==(float)-1.0) return -two25/zero; /* log1p(-1)=-inf */ - else return (x-x)/(x-x); /* log1p(x<-1)=NaN */ - } - if(ax<0x31000000) { /* |x| < 2**-29 */ - math_force_eval(two25+x); /* raise inexact */ - if (ax<0x24800000) /* |x| < 2**-54 */ - { - math_check_force_underflow (x); - return x; - } - else - return x - x*x*(float)0.5; - } - if(hx>0||hx<=((int32_t)0xbe95f61f)) { - k=0;f=x;hu=1;} /* -0.2929<x<0.41422 */ + double z = x; + uint32_t ux = asuint (x); + uint32_t ax = ux & (~0u >> 1); + if (__glibc_likely (ax < 0x3c880000)) + { + if (__glibc_unlikely (ax < 0x33000000)) + { + if (!ax) + return x; + return fmaf (x, -x, x); } - if (hx >= 0x7f800000) return x+x; - if(k!=0) { - if(hx<0x5a000000) { - u = (float)1.0+x; - GET_FLOAT_WORD(hu,u); - k = (hu>>23)-127; - /* correction term */ - c = (k>0)? (float)1.0-(u-x):x-(u-(float)1.0); - c /= u; - } else { - u = x; - GET_FLOAT_WORD(hu,u); - k = (hu>>23)-127; - c = 0; - } - hu &= 0x007fffff; - if(hu<0x3504f7) { - SET_FLOAT_WORD(u,hu|0x3f800000);/* normalize u */ - } else { - k += 1; - SET_FLOAT_WORD(u,hu|0x3f000000); /* normalize u/2 */ - hu = (0x00800000-hu)>>2; + double z2 = z * z, z4 = z2 * z2; + double f = z2 + * ((b[1] + z * b[2]) + z2 * (b[3] + z * b[4]) + + z4 * ((b[5] + z * b[6]) + z2 * b[7])); + double r = z + f; + if (__glibc_unlikely ((asuint64 (r) & 0xfffffffll) == 0)) + r += 0x1p14 * (f + (z - r)); + return r; + } + else + { + if (__glibc_unlikely (ux >= 0xbf800000u || ax >= 0x7f800000)) + return as_special (x); + uint64_t tp = asuint64 (z + 1); + int e = tp >> 52; + uint64_t m52 = tp & (~(uint64_t) 0 >> 12); + unsigned int j = (tp >> (52 - 5)) & 31; + e -= 0x3ff; + double xd = asdouble (m52 | ((uint64_t) 0x3ff << 52)); + z = xd * x0[j] - 1; + static const double c[] = + { + -0x1.3902c33434e7fp-43, 0x1.ffffffe1cbed5p-1, -0x1.ffffff7d1b014p-2, + 0x1.5564e0ed3613ap-2, -0x1.0012232a00d4ap-2 + }; + const double ln2 = 0x1.62e42fefa39efp-1; + double z2 = z * z, + r = (ln2 * e + lixb[j]) + + z * ((c[1] + z * c[2]) + z2 * (c[3] + z * c[4])); + float ub = r; + float lb = r + 2.2e-11; + if (__glibc_unlikely (ub != lb)) + { + double z4 = z2 * z2, + f = z + * ((b[0] + z * b[1]) + z2 * (b[2] + z * b[3]) + + z4 * ((b[4] + z * b[5]) + z2 * (b[6] + z * b[7]))); + const double ln2l = 0x1.7f7d1cf79abcap-20, ln2h = 0x1.62e4p-1; + double Lh = ln2h * e; + double Ll = ln2l * e; + double rl = f + Ll + lix[j]; + double tr = rl + Lh; + if (__glibc_unlikely ((asuint64 (tr) & 0xfffffffll) == 0)) + { + if (x == -0x1.247ab0p-6) + return -0x1.271f0ep-6f - 0x1p-31f; + if (x == -0x1.3a415ep-5) + return -0x1.407112p-5f + 0x1p-30f; + if (x == 0x1.fb035ap-2) + return 0x1.9bddc2p-2f + 0x1p-27f; + tr += 64 * (rl + (Lh - tr)); } - f = u-(float)1.0; - } - hfsq=(float)0.5*f*f; - if(hu==0) { /* |f| < 2**-20 */ - if(f==zero) { - if(k==0) return zero; - else {c += k*ln2_lo; return k*ln2_hi+c;} + else if (rl + (Lh - tr) == 0.0) + { + if (x == 0x1.b7fd86p-4) + return 0x1.a1ece2p-4f + 0x1p-29f; + if (x == -0x1.3a415ep-5) + return -0x1.407112p-5f + 0x1p-30f; + if (x == 0x1.43c7e2p-6) + return 0x1.409f80p-6f + 0x1p-31f; } - R = hfsq*(1.0f-0.66666666666666666f*f); - if(k==0) return f-R; else - return k*ln2_hi-((R-(k*ln2_lo+c))-f); + ub = tr; } - s = f/((float)2.0+f); - z = s*s; - R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7)))))); - if (k == 0) - return f - (hfsq - s * (hfsq + R)); - else - { - /* With GCC 7 when compiling with -Os the compiler warns - that c might be used uninitialized. This can't be true - because k must be 0 for c to be uninitialized and we - handled that computation earlier without using c. */ - DIAG_PUSH_NEEDS_COMMENT; - DIAG_IGNORE_Os_NEEDS_COMMENT (7, "-Wmaybe-uninitialized"); - return k * ln2_hi - ((hfsq - (s * (hfsq + R) - + (k * ln2_lo + c))) - f); - DIAG_POP_NEEDS_COMMENT; - } + return ub; + } } +libm_alias_float (__log1p, log1p) +strong_alias (__log1pf, __logp1f) +libm_alias_float (__logp1, logp1) diff --git a/sysdeps/ieee754/flt-32/s_log2p1f.c b/sysdeps/ieee754/flt-32/s_log2p1f.c new file mode 100644 index 0000000000..09e77dc08a --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_log2p1f.c @@ -0,0 +1,248 @@ +/* Correctly-rounded biased argument natural logarithm function for binary32 + value. + +Copyright (c) 2022-2024 Alexei Sibidanov. + +This file is part of the CORE-MATH project +project (file src/binary32/log2p1/log2p1f.c revision bc385c2). + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ + +#include <errno.h> +#include <math.h> +#include <math-underflow.h> +#include <libm-alias-float.h> +#include "math_config.h" + +float +__log2p1f (float x) +{ + static const double ix[] = + { + 0x1p+0, 0x1.fc07f01fcp-1, 0x1.f81f81f82p-1, + 0x1.f44659e4ap-1, 0x1.f07c1f07cp-1, 0x1.ecc07b302p-1, + 0x1.e9131abfp-1, 0x1.e573ac902p-1, 0x1.e1e1e1e1ep-1, + 0x1.de5d6e3f8p-1, 0x1.dae6076bap-1, 0x1.d77b654b8p-1, + 0x1.d41d41d42p-1, 0x1.d0cb58f6ep-1, 0x1.cd8568904p-1, + 0x1.ca4b3055ep-1, 0x1.c71c71c72p-1, 0x1.c3f8f01c4p-1, + 0x1.c0e070382p-1, 0x1.bdd2b8994p-1, 0x1.bacf914c2p-1, + 0x1.b7d6c3ddap-1, 0x1.b4e81b4e8p-1, 0x1.b2036406cp-1, + 0x1.af286bca2p-1, 0x1.ac5701ac6p-1, 0x1.a98ef606ap-1, + 0x1.a6d01a6dp-1, 0x1.a41a41a42p-1, 0x1.a16d3f97ap-1, + 0x1.9ec8e951p-1, 0x1.9c2d14ee4p-1, 0x1.99999999ap-1, + 0x1.970e4f80cp-1, 0x1.948b0fcd6p-1, 0x1.920fb49dp-1, + 0x1.8f9c18f9cp-1, 0x1.8d3018d3p-1, 0x1.8acb90f6cp-1, + 0x1.886e5f0acp-1, 0x1.861861862p-1, 0x1.83c977ab2p-1, + 0x1.818181818p-1, 0x1.7f405fd02p-1, 0x1.7d05f417ep-1, + 0x1.7ad2208ep-1, 0x1.78a4c8178p-1, 0x1.767dce434p-1, + 0x1.745d1745ep-1, 0x1.724287f46p-1, 0x1.702e05c0cp-1, + 0x1.6e1f76b44p-1, 0x1.6c16c16c2p-1, 0x1.6a13cd154p-1, + 0x1.681681682p-1, 0x1.661ec6a52p-1, 0x1.642c8590cp-1, + 0x1.623fa7702p-1, 0x1.605816058p-1, 0x1.5e75bb8dp-1, + 0x1.5c9882b94p-1, 0x1.5ac056b02p-1, 0x1.58ed23082p-1, + 0x1.571ed3c5p-1, 0x1.555555556p-1, 0x1.5390948f4p-1, + 0x1.51d07eae2p-1, 0x1.501501502p-1, 0x1.4e5e0a73p-1, + 0x1.4cab88726p-1, 0x1.4afd6a052p-1, 0x1.49539e3b2p-1, + 0x1.47ae147aep-1, 0x1.460cbc7f6p-1, 0x1.446f86562p-1, + 0x1.42d6625d6p-1, 0x1.414141414p-1, 0x1.3fb013fbp-1, + 0x1.3e22cbce4p-1, 0x1.3c995a47cp-1, 0x1.3b13b13b2p-1, + 0x1.3991c2c18p-1, 0x1.381381382p-1, 0x1.3698df3dep-1, + 0x1.3521cfb2cp-1, 0x1.33ae45b58p-1, 0x1.323e34a2cp-1, + 0x1.30d19013p-1, 0x1.2f684bda2p-1, 0x1.2e025c04cp-1, + 0x1.2c9fb4d82p-1, 0x1.2b404ad02p-1, 0x1.29e4129e4p-1, + 0x1.288b01288p-1, 0x1.27350b882p-1, 0x1.25e22708p-1, + 0x1.24924924ap-1, 0x1.23456789ap-1, 0x1.21fb78122p-1, + 0x1.20b470c68p-1, 0x1.1f7047dc2p-1, 0x1.1e2ef3b4p-1, + 0x1.1cf06ada2p-1, 0x1.1bb4a4046p-1, 0x1.1a7b9611ap-1, + 0x1.19453808cp-1, 0x1.181181182p-1, 0x1.16e068942p-1, + 0x1.15b1e5f76p-1, 0x1.1485f0e0ap-1, 0x1.135c81136p-1, + 0x1.12358e75ep-1, 0x1.111111112p-1, 0x1.0fef010fep-1, + 0x1.0ecf56be6p-1, 0x1.0db20a89p-1, 0x1.0c9714fbcp-1, + 0x1.0b7e6ec26p-1, 0x1.0a6810a68p-1, 0x1.0953f3902p-1, + 0x1.084210842p-1, 0x1.073260a48p-1, 0x1.0624dd2f2p-1, + 0x1.05197f7d8p-1, 0x1.041041042p-1, 0x1.03091b52p-1, + 0x1.020408102p-1, 0x1.01010101p-1, 0x1p-1 + }; + + static const double lix[] = { + 0x0p+0, -0x1.6fe50b6f1eafap-7, -0x1.6e79685c160d5p-6, + -0x1.11cd1d51955bap-5, -0x1.6bad37591e03p-5, -0x1.c4dfab908ddb5p-5, + -0x1.0eb389fab4795p-4, -0x1.3aa2fdd26ae99p-4, -0x1.663f6faca846bp-4, + -0x1.918a16e4cb157p-4, -0x1.bc84240a78a13p-4, -0x1.e72ec1181cfb1p-4, + -0x1.08c588cd964e4p-3, -0x1.1dcd19759f2e3p-3, -0x1.32ae9e27627c6p-3, + -0x1.476a9f989a58ap-3, -0x1.5c01a39fa6533p-3, -0x1.70742d4eed455p-3, + -0x1.84c2bd02d6434p-3, -0x1.98edd077e9f0ap-3, -0x1.acf5e2db31eeap-3, + -0x1.c0db6cddaa82dp-3, -0x1.d49ee4c33121ap-3, -0x1.e840be751d775p-3, + -0x1.fbc16b9003e0bp-3, -0x1.0790adbae3fcp-2, -0x1.11307dad465b5p-2, + -0x1.1ac05b2924cc5p-2, -0x1.24407ab0cc41p-2, -0x1.2db10fc4ea424p-2, + -0x1.37124cea58697p-2, -0x1.406463b1d455dp-2, -0x1.49a784bcbaa37p-2, + -0x1.52dbdfc4f341dp-2, -0x1.5c01a39ff2c9bp-2, -0x1.6518fe46abaa5p-2, + -0x1.6e221cd9d6933p-2, -0x1.771d2ba7f5791p-2, -0x1.800a56315ee2ap-2, + -0x1.88e9c72df8611p-2, -0x1.91bba891d495fp-2, -0x1.9a8023920fa4dp-2, + -0x1.a33760a7fbca6p-2, -0x1.abe18797d2effp-2, -0x1.b47ebf734b923p-2, + -0x1.bd0f2e9eb2b84p-2, -0x1.c592fad2be1aap-2, -0x1.ce0a4923cf5e6p-2, + -0x1.d6753e02f4ebcp-2, -0x1.ded3fd445afp-2, -0x1.e726aa1e558fep-2, + -0x1.ef6d67325ba38p-2, -0x1.f7a8568c8aea6p-2, -0x1.ffd799a81be87p-2, + 0x1.f804ae8d33c4p-2, 0x1.efec61b04af4ep-2, 0x1.e7df5fe572606p-2, + 0x1.dfdd89d5b0009p-2, 0x1.d7e6c0abbd924p-2, 0x1.cffae611a74d6p-2, + 0x1.c819dc2d8578cp-2, 0x1.c043859e5bdbcp-2, 0x1.b877c57b47c04p-2, + 0x1.b0b67f4f29a66p-2, 0x1.a8ff97183ed07p-2, 0x1.a152f14293c74p-2, + 0x1.99b072a9289cap-2, 0x1.921800927e284p-2, 0x1.8a8980ac4113p-2, + 0x1.8304d90c2859dp-2, 0x1.7b89f02cbd49ap-2, 0x1.7418aceb84ab1p-2, + 0x1.6cb0f68656c95p-2, 0x1.6552b49993dc2p-2, 0x1.5dfdcf1eacd7bp-2, + 0x1.56b22e6b97c18p-2, 0x1.4f6fbb2ce6943p-2, 0x1.48365e6957b42p-2, + 0x1.4106017c0dbcfp-2, 0x1.39de8e15727d9p-2, 0x1.32bfee37489bcp-2, + 0x1.2baa0c34989c3p-2, 0x1.249cd2b177fd5p-2, 0x1.1d982c9d50468p-2, + 0x1.169c0536677acp-2, 0x1.0fa848045f67bp-2, 0x1.08bce0d9a7c6p-2, + 0x1.01d9bbcf66a2cp-2, 0x1.f5fd8a90e2d85p-3, 0x1.e857d3d3af1e5p-3, + 0x1.dac22d3ec5f4ep-3, 0x1.cd3c712db459ap-3, 0x1.bfc67a7ff3c22p-3, + 0x1.b2602497678f4p-3, 0x1.a5094b555a1f8p-3, 0x1.97c1cb136b96fp-3, + 0x1.8a8980ac8652dp-3, 0x1.7d60496c83f66p-3, 0x1.7046031c7cdafp-3, + 0x1.633a8bf460335p-3, 0x1.563dc2a08b102p-3, 0x1.494f863bbc1dep-3, + 0x1.3c6fb6507a37ep-3, 0x1.2f9e32d5257ecp-3, 0x1.22dadc2a627efp-3, + 0x1.1625931802e49p-3, 0x1.097e38cef9519p-3, 0x1.f9c95dc138295p-4, + 0x1.e0b1ae90505f6p-4, 0x1.c7b528b5fcffap-4, 0x1.aed391abb17a1p-4, + 0x1.960caf9bd35eap-4, 0x1.7d60496e3edebp-4, 0x1.64ce26bf2108ep-4, + 0x1.4c560fe5b573bp-4, 0x1.33f7cde24adfbp-4, 0x1.1bb32a5ed9353p-4, + 0x1.0387efbd3006ep-4, 0x1.d6ebd1f1d0955p-5, 0x1.a6f9c37a8beabp-5, + 0x1.77394c9d6762cp-5, 0x1.47aa07358e1a4p-5, 0x1.184b8e4d490efp-5, + 0x1.d23afc4d95c78p-6, 0x1.743ee8678a7cbp-6, 0x1.16a21e243bf78p-6, + 0x1.72c7ba20c907ep-7, 0x1.720d9c0536e17p-8, 0x0p+0 + }; + + double z = x; + uint32_t ux = asuint (x); + uint32_t ax = ux & (~0u >> 1); + if (__glibc_unlikely (ux >= 0x17fu << 23)) + { /* x <= -1 */ + if (ux == (0x17fu << 23)) + return __math_divzerof (1); + if (ux > (0x1ffu << 23)) + return x + x; /* nan */ + return __math_invalidf (x); + } + else if (__glibc_unlikely (ax >= (0xff << 23))) + { /* +inf, nan */ + if (ax > (0xff << 23)) + return x + x; /* nan */ + return INFINITY; + } + else if (__glibc_likely (ax < 0x3cb7aa26u)) + { /* |x| < 0x1.6f544cp-6 */ + double z2 = z * z, z4 = z2 * z2; + if ( __glibc_likely (ax < 0x3b9d9d34u)) + { /* |x| < 0x1.3b3a68p-8 */ + if (__glibc_likely (ax < 0x39638a7eu)) + { /* |x| < 0x1.c714fcp-13 */ + if (__glibc_likely (ax < 0x329c5639u)) + { /* |x| < 0x1.38ac72p-26 */ + static const double c[] = + { + 0x1.71547652b82fep+0, -0x1.71547652b82ffp-1 + }; + return z * (c[0] + z * c[1]); + } + else + { + if (__glibc_unlikely (ux == 0x32ff7045u)) + return 0x1.70851ap-25f - 0x1.8p-80f; + if (__glibc_unlikely (ux == 0xb395efbbu)) + return -0x1.b0a00ap-24f + 0x1p-76f; + if (__glibc_unlikely (ux == 0x35a14df7u)) + return 0x1.d16d2p-20f + 0x1p-72f; + if (__glibc_unlikely (ux == 0x3841cb81u)) + return 0x1.17949ep-14f + 0x1p-67f; + static const double c[] = + { + 0x1.71547652b82fep+0, -0x1.71547652b82fdp-1, + 0x1.ec709ead0c9a7p-2, -0x1.7154773c1cb29p-2 + }; + return z * ((c[0] + z * c[1]) + z2 * (c[2] + z * c[3])); + } + } + else + { + if (__glibc_unlikely (ux == 0xbac9363du)) + return -0x1.2282aap-9f + 0x1p-61f; + static const double c[] = + { + 0x1.71547652b82fep+0, -0x1.71547652b83p-1, + 0x1.ec709dc28f51bp-2, -0x1.7154765157748p-2, + 0x1.2778a510a3682p-2, -0x1.ec745df1551fcp-3 + }; + return z + * ((c[0] + z * c[1]) + z2 * (c[2] + z * c[3]) + + z4 * ((c[4] + z * c[5]))); + } + } + else + { + static const double c[] = + { + 0x1.71547652b82fep+0, -0x1.71547652b82fbp-1, + 0x1.ec709dc3b6a73p-2, -0x1.71547652dc09p-2, + 0x1.2776c1a88901p-2, -0x1.ec7095bd4d208p-3, + 0x1.a66bec7fc8f7p-3, -0x1.71a900fc3f3f9p-3 + }; + return z + * ((c[0] + z * c[1]) + z2 * (c[2] + z * c[3]) + + z4 * ((c[4] + z * c[5]) + z2 * (c[6] + z * c[7]))); + } + } + else + { /* |x| >= 0x1.6f544cp-6 */ + float h, l; + /* With gcc 6.3.0, if we return 0x1.e90026p+4f + 0x1.fp-21 + in the second exceptional case, with rounding up it yields + 0x1.e90026p+4 which is incorrect, thus we use this workaround. See + https://gcc.gnu.org/bugzilla/show_bug.cgi?id=112367. */ + if (__glibc_unlikely (ux == 0x52928e33u)) + { + h = 0x1.318ffap+5f; + l = 0x1.fp-20f; + return h + l; + } + if (__glibc_unlikely (ux == 0x4ebd09e3u)) + { + h = 0x1.e90026p+4f; + l = 0x1.fp-21; + return h + l; + } + uint64_t tp = asuint64 (z + 1.0); + uint64_t m = tp & (~(uint64_t) 0 >> 12); + int e = (tp >> 52) - 0x3ff; + int j = (m + ((int64_t) 1 << (52 - 8))) >> (52 - 7), k = j > 53; + e += k; + double xd = asdouble (m | (uint64_t) 0x3ff << 52); + z = fma (xd, ix[j], -1.0); + static const double c[] = + { + 0x1.71547652b82fep+0, -0x1.71547652b82ffp-1, 0x1.ec709dc32988bp-2, + -0x1.715476521ec2bp-2, 0x1.277801a1ad904p-2, -0x1.ec731704d6a88p-3 + }; + double z2 = z * z; + double c0 = c[0] + z * c[1]; + double c2 = c[2] + z * c[3]; + double c4 = c[4] + z * c[5]; + c0 += z2 * (c2 + z2 * c4); + return (z * c0 - lix[j]) + e; + } +} +libm_alias_float (__log2p1, log2p1) diff --git a/sysdeps/ieee754/flt-32/w_log1pf.c b/sysdeps/ieee754/flt-32/w_log1pf.c new file mode 100644 index 0000000000..1cc8931700 --- /dev/null +++ b/sysdeps/ieee754/flt-32/w_log1pf.c @@ -0,0 +1 @@ +/* Not needed. */ |