diff options
30 files changed, 357 insertions, 601 deletions
diff --git a/SHARED-FILES b/SHARED-FILES index e6e29bcadc..033ce7f092 100644 --- a/SHARED-FILES +++ b/SHARED-FILES @@ -280,3 +280,11 @@ sysdeps/ieee754/flt-32/s_erfcf.c (file src/binary32/erfc/erfcf.c in CORE-MATH) - The code was adapted to use glibc code style and internal functions to handle errno, overflow, and underflow. +sysdeps/ieee754/flt-32/e_lgammaf_r.c: + (file src/binary32/lgamma/lgammaf.c in CORE-MATH) + - change the function name from cr_lgammaf to __ieee754_lgammaf_r + - add "int *signgamp" as 2nd argument and add at the beginning: + if (signgamp != NULL) *signgamp = 1; + - remove the errno stuff (this is done by the wrapper) + - replace 0x1p127f * 0x1p127f by math_narrow_eval (x * 0x1p127f) + - add libm_alias_finite (__ieee754_lgammaf_r, __lgammaf_r) at the end diff --git a/sysdeps/aarch64/libm-test-ulps b/sysdeps/aarch64/libm-test-ulps index fc10f7f80d..1d3d1f9b6a 100644 --- a/sysdeps/aarch64/libm-test-ulps +++ b/sysdeps/aarch64/libm-test-ulps @@ -1288,22 +1288,18 @@ ldouble: 7 Function: "lgamma": double: 3 -float: 4 ldouble: 5 Function: "lgamma_downward": double: 4 -float: 4 ldouble: 8 Function: "lgamma_towardzero": double: 4 -float: 3 ldouble: 5 Function: "lgamma_upward": double: 4 -float: 5 ldouble: 8 Function: "log": diff --git a/sysdeps/alpha/fpu/libm-test-ulps b/sysdeps/alpha/fpu/libm-test-ulps index 3678bc38e3..7256e674bb 100644 --- a/sysdeps/alpha/fpu/libm-test-ulps +++ b/sysdeps/alpha/fpu/libm-test-ulps @@ -1134,22 +1134,18 @@ ldouble: 7 Function: "lgamma": double: 4 -float: 7 ldouble: 5 Function: "lgamma_downward": double: 5 -float: 7 ldouble: 8 Function: "lgamma_towardzero": double: 5 -float: 6 ldouble: 5 Function: "lgamma_upward": double: 5 -float: 6 ldouble: 8 Function: "log": diff --git a/sysdeps/arc/fpu/libm-test-ulps b/sysdeps/arc/fpu/libm-test-ulps index 68d67b6a67..66a2b541c6 100644 --- a/sysdeps/arc/fpu/libm-test-ulps +++ b/sysdeps/arc/fpu/libm-test-ulps @@ -917,19 +917,15 @@ float: 9 Function: "lgamma": double: 7 -float: 6 Function: "lgamma_downward": double: 6 -float: 5 Function: "lgamma_towardzero": double: 7 -float: 6 Function: "lgamma_upward": double: 7 -float: 6 Function: "log": double: 1 diff --git a/sysdeps/arc/nofpu/libm-test-ulps b/sysdeps/arc/nofpu/libm-test-ulps index dc2499b56a..38836ddc38 100644 --- a/sysdeps/arc/nofpu/libm-test-ulps +++ b/sysdeps/arc/nofpu/libm-test-ulps @@ -223,7 +223,6 @@ float: 4 Function: "lgamma": double: 4 -float: 7 Function: "log10": double: 2 diff --git a/sysdeps/arm/libm-test-ulps b/sysdeps/arm/libm-test-ulps index 000a5af492..2651046cfa 100644 --- a/sysdeps/arm/libm-test-ulps +++ b/sysdeps/arm/libm-test-ulps @@ -911,19 +911,15 @@ float: 5 Function: "lgamma": double: 4 -float: 7 Function: "lgamma_downward": double: 5 -float: 7 Function: "lgamma_towardzero": double: 5 -float: 6 Function: "lgamma_upward": double: 5 -float: 6 Function: "log": float: 1 diff --git a/sysdeps/csky/fpu/libm-test-ulps b/sysdeps/csky/fpu/libm-test-ulps index ed373cd353..02b4cb4934 100644 --- a/sysdeps/csky/fpu/libm-test-ulps +++ b/sysdeps/csky/fpu/libm-test-ulps @@ -875,19 +875,15 @@ float: 5 Function: "lgamma": double: 4 -float: 7 Function: "lgamma_downward": double: 5 -float: 7 Function: "lgamma_towardzero": double: 5 -float: 6 Function: "lgamma_upward": double: 5 -float: 6 Function: "log10": double: 2 diff --git a/sysdeps/csky/nofpu/libm-test-ulps b/sysdeps/csky/nofpu/libm-test-ulps index 48a8c7351f..34312f5a06 100644 --- a/sysdeps/csky/nofpu/libm-test-ulps +++ b/sysdeps/csky/nofpu/libm-test-ulps @@ -873,19 +873,15 @@ float: 5 Function: "lgamma": double: 4 -float: 7 Function: "lgamma_downward": double: 5 -float: 4 Function: "lgamma_towardzero": double: 5 -float: 4 Function: "lgamma_upward": double: 5 -float: 5 Function: "log": float: 1 diff --git a/sysdeps/hppa/fpu/libm-test-ulps b/sysdeps/hppa/fpu/libm-test-ulps index 9c298ac5e3..976661541b 100644 --- a/sysdeps/hppa/fpu/libm-test-ulps +++ b/sysdeps/hppa/fpu/libm-test-ulps @@ -934,20 +934,16 @@ float: 5 Function: "lgamma": double: 4 -float: 7 ldouble: 1 Function: "lgamma_downward": double: 5 -float: 7 Function: "lgamma_towardzero": double: 5 -float: 6 Function: "lgamma_upward": double: 5 -float: 6 Function: "log": double: 1 diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps index 5284078bee..170e7cfc65 100644 --- a/sysdeps/i386/fpu/libm-test-ulps +++ b/sysdeps/i386/fpu/libm-test-ulps @@ -1354,25 +1354,21 @@ ldouble: 5 Function: "lgamma": double: 4 -float: 5 float128: 5 ldouble: 4 Function: "lgamma_downward": double: 5 -float: 5 float128: 8 ldouble: 7 Function: "lgamma_towardzero": double: 5 -float: 6 float128: 5 ldouble: 7 Function: "lgamma_upward": double: 5 -float: 6 float128: 8 ldouble: 6 diff --git a/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps b/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps index da00d80ba7..a9ce632e6a 100644 --- a/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps +++ b/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps @@ -1357,7 +1357,6 @@ ldouble: 5 Function: "lgamma": double: 4 -float: 5 float128: 5 ldouble: 4 diff --git a/sysdeps/ieee754/flt-32/e_lgammaf_r.c b/sysdeps/ieee754/flt-32/e_lgammaf_r.c index a1a3a60454..cb65513056 100644 --- a/sysdeps/ieee754/flt-32/e_lgammaf_r.c +++ b/sysdeps/ieee754/flt-32/e_lgammaf_r.c @@ -1,247 +1,367 @@ -/* e_lgammaf_r.c -- float version of e_lgamma_r.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 logarithm of the absolute value of the gamma function + for binary32 value. +Copyright (c) 2023, 2024 Alexei Sibidanov. + +This file is part of the CORE-MATH project +project (file src/binary32/lgamma/lgammaf.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. +*/ + +/* Changes with respect to the original CORE-MATH code: + - removed the dealing with errno + (this is done in the wrapper math/w_lgammaf_compat2.c). + - usage of math_narrow_eval to deal with underflow/overflow. + - deal with signamp. */ + +#include <array_length.h> +#include <stdint.h> #include <math.h> -#include <math-narrow-eval.h> -#include <math_private.h> -#include <libc-diag.h> #include <libm-alias-finite.h> +#include <limits.h> +#include <math-narrow-eval.h> +#include "math_config.h" -static const float -two23= 8.3886080000e+06, /* 0x4b000000 */ -half= 5.0000000000e-01, /* 0x3f000000 */ -one = 1.0000000000e+00, /* 0x3f800000 */ -pi = 3.1415927410e+00, /* 0x40490fdb */ -a0 = 7.7215664089e-02, /* 0x3d9e233f */ -a1 = 3.2246702909e-01, /* 0x3ea51a66 */ -a2 = 6.7352302372e-02, /* 0x3d89f001 */ -a3 = 2.0580807701e-02, /* 0x3ca89915 */ -a4 = 7.3855509982e-03, /* 0x3bf2027e */ -a5 = 2.8905137442e-03, /* 0x3b3d6ec6 */ -a6 = 1.1927076848e-03, /* 0x3a9c54a1 */ -a7 = 5.1006977446e-04, /* 0x3a05b634 */ -a8 = 2.2086278477e-04, /* 0x39679767 */ -a9 = 1.0801156895e-04, /* 0x38e28445 */ -a10 = 2.5214456400e-05, /* 0x37d383a2 */ -a11 = 4.4864096708e-05, /* 0x383c2c75 */ -tc = 1.4616321325e+00, /* 0x3fbb16c3 */ -tf = -1.2148628384e-01, /* 0xbdf8cdcd */ -/* tt = -(tail of tf) */ -tt = 6.6971006518e-09, /* 0x31e61c52 */ -t0 = 4.8383611441e-01, /* 0x3ef7b95e */ -t1 = -1.4758771658e-01, /* 0xbe17213c */ -t2 = 6.4624942839e-02, /* 0x3d845a15 */ -t3 = -3.2788541168e-02, /* 0xbd064d47 */ -t4 = 1.7970675603e-02, /* 0x3c93373d */ -t5 = -1.0314224288e-02, /* 0xbc28fcfe */ -t6 = 6.1005386524e-03, /* 0x3bc7e707 */ -t7 = -3.6845202558e-03, /* 0xbb7177fe */ -t8 = 2.2596477065e-03, /* 0x3b141699 */ -t9 = -1.4034647029e-03, /* 0xbab7f476 */ -t10 = 8.8108185446e-04, /* 0x3a66f867 */ -t11 = -5.3859531181e-04, /* 0xba0d3085 */ -t12 = 3.1563205994e-04, /* 0x39a57b6b */ -t13 = -3.1275415677e-04, /* 0xb9a3f927 */ -t14 = 3.3552918467e-04, /* 0x39afe9f7 */ -u0 = -7.7215664089e-02, /* 0xbd9e233f */ -u1 = 6.3282704353e-01, /* 0x3f2200f4 */ -u2 = 1.4549225569e+00, /* 0x3fba3ae7 */ -u3 = 9.7771751881e-01, /* 0x3f7a4bb2 */ -u4 = 2.2896373272e-01, /* 0x3e6a7578 */ -u5 = 1.3381091878e-02, /* 0x3c5b3c5e */ -v1 = 2.4559779167e+00, /* 0x401d2ebe */ -v2 = 2.1284897327e+00, /* 0x4008392d */ -v3 = 7.6928514242e-01, /* 0x3f44efdf */ -v4 = 1.0422264785e-01, /* 0x3dd572af */ -v5 = 3.2170924824e-03, /* 0x3b52d5db */ -s0 = -7.7215664089e-02, /* 0xbd9e233f */ -s1 = 2.1498242021e-01, /* 0x3e5c245a */ -s2 = 3.2577878237e-01, /* 0x3ea6cc7a */ -s3 = 1.4635047317e-01, /* 0x3e15dce6 */ -s4 = 2.6642270386e-02, /* 0x3cda40e4 */ -s5 = 1.8402845599e-03, /* 0x3af135b4 */ -s6 = 3.1947532989e-05, /* 0x3805ff67 */ -r1 = 1.3920053244e+00, /* 0x3fb22d3b */ -r2 = 7.2193557024e-01, /* 0x3f38d0c5 */ -r3 = 1.7193385959e-01, /* 0x3e300f6e */ -r4 = 1.8645919859e-02, /* 0x3c98bf54 */ -r5 = 7.7794247773e-04, /* 0x3a4beed6 */ -r6 = 7.3266842264e-06, /* 0x36f5d7bd */ -w0 = 4.1893854737e-01, /* 0x3ed67f1d */ -w1 = 8.3333335817e-02, /* 0x3daaaaab */ -w2 = -2.7777778450e-03, /* 0xbb360b61 */ -w3 = 7.9365057172e-04, /* 0x3a500cfd */ -w4 = -5.9518753551e-04, /* 0xba1c065c */ -w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */ -w6 = -1.6309292987e-03; /* 0xbad5c4e8 */ - -static const float zero= 0.0000000000e+00; - -static float -sin_pif(float x) +static double +as_r7 (double x, const double *c) { - float y,z; - int n,ix; - - GET_FLOAT_WORD(ix,x); - ix &= 0x7fffffff; - - if(ix<0x3e800000) return __sinf (pi*x); - y = -x; /* x is assume negative */ - - /* - * argument reduction, make sure inexact flag not raised if input - * is an integer - */ - z = floorf(y); - if(z!=y) { /* inexact anyway */ - y *= (float)0.5; - y = (float)2.0*(y - floorf(y)); /* y = |x| mod 2.0 */ - n = (int) (y*(float)4.0); - } else { - if(ix>=0x4b800000) { - y = zero; n = 0; /* y must be even */ - } else { - if(ix<0x4b000000) z = y+two23; /* exact */ - GET_FLOAT_WORD(n,z); - n &= 1; - y = n; - n<<= 2; - } - } - switch (n) { - case 0: y = __sinf (pi*y); break; - case 1: - case 2: y = __cosf (pi*((float)0.5-y)); break; - case 3: - case 4: y = __sinf (pi*(one-y)); break; - case 5: - case 6: y = -__cosf (pi*(y-(float)1.5)); break; - default: y = __sinf (pi*(y-(float)2.0)); break; - } - return -y; + return (((x - c[0]) * (x - c[1])) * ((x - c[2]) * (x - c[3]))) + * (((x - c[4]) * (x - c[5])) * ((x - c[6]))); } +static double +as_r8 (double x, const double *c) +{ + return (((x - c[0]) * (x - c[1])) * ((x - c[2]) * (x - c[3]))) + * (((x - c[4]) * (x - c[5])) * ((x - c[6]) * (x - c[7]))); +} + +static double +as_sinpi (double x) +{ + static const double c[] = + { + 0x1p+2, -0x1.de9e64df22ea4p+1, 0x1.472be122401f8p+0, + -0x1.d4fcd82df91bp-3, 0x1.9f05c97e0aab2p-6, -0x1.f3091c427b611p-10, + 0x1.b22c9bfdca547p-14, -0x1.15484325ef569p-18 + }; + x -= 0.5; + double x2 = x * x, x4 = x2 * x2, x8 = x4 * x4; + return (0.25 - x2) + * ((c[0] + x2 * c[1]) + x4 * (c[2] + x2 * c[3]) + + x8 * ((c[4] + x2 * c[5]) + x4 * (c[6] + x2 * c[7]))); +} + +static double +as_ln (double x) +{ + uint64_t t = asuint64 (x); + int e = (t >> 52) - 0x3ff; + static const double c[] = + { + 0x1.fffffffffff24p-1, -0x1.ffffffffd1d67p-2, 0x1.55555537802dep-2, + -0x1.ffffeca81b866p-3, 0x1.999611761d772p-3, -0x1.54f3e581b61bfp-3, + 0x1.1e642b4cb5143p-3, -0x1.9115a5af1e1edp-4 + }; + static const double il[] = + { + 0x1.59caeec280116p-57, 0x1.f0a30c01162aap-5, 0x1.e27076e2af2ebp-4, + 0x1.5ff3070a793d6p-3, 0x1.c8ff7c79a9a2p-3, 0x1.1675cababa60fp-2, + 0x1.4618bc21c5ec2p-2, 0x1.739d7f6bbd007p-2, 0x1.9f323ecbf984dp-2, + 0x1.c8ff7c79a9a21p-2, 0x1.f128f5faf06ecp-2, 0x1.0be72e4252a83p-1, + 0x1.1e85f5e7040d1p-1, 0x1.307d7334f10bep-1, 0x1.41d8fe84672afp-1, + 0x1.52a2d265bc5abp-1 + }; + static const double ix[] = + { + 0x1p+0, 0x1.e1e1e1e1e1e1ep-1, 0x1.c71c71c71c71cp-1, + 0x1.af286bca1af28p-1, 0x1.999999999999ap-1, 0x1.8618618618618p-1, + 0x1.745d1745d1746p-1, 0x1.642c8590b2164p-1, 0x1.5555555555555p-1, + 0x1.47ae147ae147bp-1, 0x1.3b13b13b13b14p-1, 0x1.2f684bda12f68p-1, + 0x1.2492492492492p-1, 0x1.1a7b9611a7b96p-1, 0x1.1111111111111p-1, + 0x1.0842108421084p-1 + }; + int i = (t >> 48) & 0xf; + t = (t & (~UINT64_C(0) >> 12)) | (INT64_C(0x3ff) << 52); + double z = ix[i] * asdouble (t) - 1; + double z2 = z * z, z4 = z2 * z2; + return e * 0x1.62e42fefa39efp-1 + il[i] + + z * ((c[0] + z * c[1]) + z2 * (c[2] + z * c[3]) + + z4 * ((c[4] + z * c[5]) + z2 * (c[6] + z * c[7]))); +} float -__ieee754_lgammaf_r(float x, int *signgamp) +__ieee754_lgammaf_r (float x, int *signgamp) { - float t,y,z,nadj,p,p1,p2,p3,q,r,w; - int i,hx,ix; - - GET_FLOAT_WORD(hx,x); - - /* purge off +-inf, NaN, +-0, and negative arguments */ - *signgamp = 1; - ix = hx&0x7fffffff; - if(__builtin_expect(ix>=0x7f800000, 0)) return x*x; - if(__builtin_expect(ix==0, 0)) - { - if (hx < 0) - *signgamp = -1; - return one/fabsf(x); - } - if(__builtin_expect(ix<0x30800000, 0)) { - /* |x|<2**-30, return -log(|x|) */ - if(hx<0) { - *signgamp = -1; - return -__ieee754_logf(-x); - } else return -__ieee754_logf(x); + static const struct + { + float x; + float f; + float df; + } tb[] = { + { -0x1.efc2a2p+14, -0x1.222dbcp+18, -0x1p-7 }, + { -0x1.627346p+7, -0x1.73235ep+9, -0x1p-16 }, + { -0x1.08b14p+4, -0x1.f0cbe6p+4, -0x1p-21 }, + { -0x1.69d628p+3, -0x1.0eac2ap+4, -0x1p-21 }, + { -0x1.904902p+2, -0x1.65532cp+2, 0x1p-23 }, + { -0x1.9272d2p+1, -0x1.170b98p- |