From: Adhemerval Zanella Date: Tue, 27 Jan 2026 19:56:41 +0000 (+0000) Subject: math: Sync log1pf with CORE-MATH X-Git-Url: http://git.ipfire.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3069c7951d83a3d5d4ef87ebcd54adb64257f4e5;p=thirdparty%2Fglibc.git math: Sync log1pf with CORE-MATH The performance is similar, with some minor regression: latency master sync improvement x86_64 38.2841 38.1375 0.38% x86_64v2 37.7338 37.4292 0.81% x86_64v3 31.3500 32.3576 -3.21% aarch64 13.7384 13.9030 -1.20% armhf-vpfv4 15.5730 16.5105 -6.02% powerpc64le 7.6038 7.5757 0.37% reciprocal-throughput master sync improvement x86_64 12.4910 11.9683 4.18% x86_64v2 12.2935 11.7614 4.33% x86_64v3 11.5444 10.6369 7.86% aarch64 7.7262 7.8954 -2.19% armhf-vpfv4 8.3502 8.8741 -6.27% powerpc64le 3.5883 3.5259 1.74% x86_64 / i686 gcc version 15.2.1 20260112. Ryzen 5900X aarch64: gcc version 15.2.1 20251105, Neoverse-N1 armv7a-vpfv4: gcc version 15.2.1 20251105, Neoverse-N1 powerpc64le: gcc version 14.2.1 20241230, POWER10 The sync also improves the internal table size, the s_log1pf.os 'size' output shows: size master sync improvement x86_64 2078 1641 21.03% x86_64v2 2078 1641 21.03% x86_64v3 1975 1514 23.34% aarch64 1808 1336 26.11% armhf-vpfv4 1716 1284 25.17% powerpc64le 2132 1616 24.20% Checked on aarch64-linux-gnu, arm-linux-gnueabihf, powerpc64le-linux-gnu, i686-linux-gnu, and x86_64-linux-gnu. Reviewed-by: Paul Zimmermann --- diff --git a/SHARED-FILES b/SHARED-FILES index e942ae388e..7550bc8b0b 100644 --- a/SHARED-FILES +++ b/SHARED-FILES @@ -296,7 +296,7 @@ core-math: sysdeps/ieee754/flt-32/s_expm1f.c # src/binary32/log10p1/log10p1f.c revision bc385c2 sysdeps/ieee754/flt-32/s_log10p1f.c - # src/binary32/log1p/log1pf.c revision bc385c2 + # src/binary32/log1p/log1pf.c revision 24ef43a1 sysdeps/ieee754/flt-32/s_log1pf.c # src/binary32/log2p1/log2p1f.c revision bc385c2 sysdeps/ieee754/flt-32/s_log2p1f.c diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in index 2001baa605..4fd72b3ab6 100644 --- a/math/auto-libm-test-in +++ b/math/auto-libm-test-in @@ -7650,6 +7650,7 @@ log1p -0x4.f37d3c9ce0b14bdd86eb157df5d4p-4 log1p 0x7.2eca50c4d93196362b4f37f6e8dcp-4 log1p -0x6.3fef3067427e43dfcde9e48f74bcp-4 log1p 0x6.af53d00fd2845d4772260ef5adc4p-4 +log1p -0x1.fffffcp-127 log2 1 log2 e diff --git a/math/auto-libm-test-out-log1p b/math/auto-libm-test-out-log1p index f7d3b35e6d..1b5326a2d1 100644 --- a/math/auto-libm-test-out-log1p +++ b/math/auto-libm-test-out-log1p @@ -2711,3 +2711,28 @@ log1p 0x6.af53d00fd2845d4772260ef5adc4p-4 = log1p tonearest ibm128 0x6.af53d00fd2845d4772260ef5acp-4 : 0x5.95f3ec4683fa14a354007a53e8p-4 : inexact-ok = log1p towardzero ibm128 0x6.af53d00fd2845d4772260ef5acp-4 : 0x5.95f3ec4683fa14a354007a53e8p-4 : inexact-ok = log1p upward ibm128 0x6.af53d00fd2845d4772260ef5acp-4 : 0x5.95f3ec4683fa14a354007a53eap-4 : inexact-ok +log1p -0x1.fffffcp-127 += log1p downward binary32 -0x3.fffff8p-128 : -0x4p-128 : inexact-ok underflow-ok errno-erange-ok += log1p tonearest binary32 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok underflow-ok errno-erange-ok += log1p towardzero binary32 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok underflow-ok errno-erange-ok += log1p upward binary32 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok underflow-ok errno-erange-ok += log1p downward binary64 -0x3.fffff8p-128 : -0x3.fffff80000002p-128 : inexact-ok += log1p tonearest binary64 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok += log1p towardzero binary64 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok += log1p upward binary64 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok += log1p downward intel96 -0x3.fffff8p-128 : -0x3.fffff80000000004p-128 : inexact-ok += log1p tonearest intel96 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok += log1p towardzero intel96 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok += log1p upward intel96 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok += log1p downward m68k96 -0x3.fffff8p-128 : -0x3.fffff80000000004p-128 : inexact-ok += log1p tonearest m68k96 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok += log1p towardzero m68k96 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok += log1p upward m68k96 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok += log1p downward binary128 -0x3.fffff8p-128 : -0x3.fffff80000000000000000000002p-128 : inexact-ok += log1p tonearest binary128 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok += log1p towardzero binary128 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok += log1p upward binary128 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok += log1p downward ibm128 -0x3.fffff8p-128 : -0x3.fffff800000000000000000001p-128 : inexact-ok += log1p tonearest ibm128 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok += log1p towardzero ibm128 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok += log1p upward ibm128 -0x3.fffff8p-128 : -0x3.fffff8p-128 : inexact-ok diff --git a/sysdeps/ieee754/flt-32/s_log1pf.c b/sysdeps/ieee754/flt-32/s_log1pf.c index f03c60df34..4cd5fa10aa 100644 --- a/sysdeps/ieee754/flt-32/s_log1pf.c +++ b/sysdeps/ieee754/flt-32/s_log1pf.c @@ -4,7 +4,7 @@ Copyright (c) 2023, 2024 Alexei Sibidanov. This file is part of the CORE-MATH project -project (file src/binary32/log1p/log1pf.c revision bc385c2). +project (file src/binary32/log1p/log1pf.c revision 24ef43a1). Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal @@ -48,45 +48,33 @@ as_special (float x) float __log1pf (float x) { + /* the reciprocal 1/(1+(j+0.5)/32) is rounded to 23 bits */ 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 + { 0x1.f81f80p-1, 0x1.e9131cp-1, 0x1.dae608p-1, 0x1.cd8568p-1, + 0x1.c0e070p-1, 0x1.b4e81cp-1, 0x1.a98ef8p-1, 0x1.9ec8e8p-1, + 0x1.948b10p-1, 0x1.8acb90p-1, 0x1.818180p-1, 0x1.78a4c8p-1, + 0x1.702e04p-1, 0x1.681680p-1, 0x1.605818p-1, 0x1.58ed24p-1, + 0x1.51d080p-1, 0x1.4afd6cp-1, 0x1.446f88p-1, 0x1.3e22ccp-1, + 0x1.381380p-1, 0x1.323e34p-1, 0x1.2c9fb4p-1, 0x1.27350cp-1, + 0x1.21fb78p-1, 0x1.1cf06cp-1, 0x1.181180p-1, 0x1.135c80p-1, + 0x1.0ecf58p-1, 0x1.0a6810p-1, 0x1.0624dcp-1, 0x1.020408p-1 }; + + /* the logarithm of the reciprocal is offset by 0x1.7654p-37 so + log1p_fast(x) - log1p(x) > 0 */ + static const double lix[] = { + 0x1.fc0b0b1599ce4p-7, 0x1.77457a64a42abp-5, 0x1.341d74627847dp-4, + 0x1.a926d8a568810p-4, 0x1.0d77e8cd667aap-3, 0x1.44d2b38d15679p-3, + 0x1.7ab886a16b2b3p-3, 0x1.af3c9b686996dp-3, 0x1.e27075e30cc37p-3, + 0x1.0a3250a767d98p-2, 0x1.229423bd2662ep-2, 0x1.3a64c596c3292p-2, + 0x1.51aadd530e505p-2, 0x1.686c85e9e0177p-2, 0x1.7eaf7df859caep-2, + 0x1.94793ee2403b3p-2, 0x1.a9cec5a9cf512p-2, 0x1.beb4d3baa086fp-2, + 0x1.d32fe2a03d8b5p-2, 0x1.e744257d97431p-2, 0x1.faf58cf7bdfe7p-2, + 0x1.0723e6d1e5599p-1, 0x1.109f3b52ec2f3p-1, 0x1.19ee6a7693fc5p-1, + 0x1.23130d9c03597p-1, 0x1.2c0e9cc4604f1p-1, 0x1.34e28bd9e5837p-1, + 0x1.3d9028a72cd5fp-1, 0x1.4618b9c1dd52dp-1, 0x1.4e7d825b8d20bp-1, + 0x1.56bf9fab56a02p-1, 0x1.5ee02ab258ccap-1 + }; static const double b[] = { 0x1p+0, @@ -98,14 +86,23 @@ __log1pf (float x) 0x1.24adeca50e2bcp-3, -0x1.001ba33bf57cfp-3 }; + static const double c[] = + { + 0x1.ffffffe1eac82p-1, -0x1.ffffff7da1724p-2, 0x1.5564d8fa59d0cp-2, + -0x1.001219d3dba2ap-2 + }; double z = x; uint32_t ux = asuint (x); + if (__glibc_unlikely (ux >= 0xbf800000u)) + return as_special (x); // x<=-1, x=-inf, x=-nan uint32_t ax = ux & (~0u >> 1); - if (__glibc_likely (ax < 0x3c880000)) - { - if (__glibc_unlikely (ax < 0x33000000)) - { + if (__glibc_unlikely (ax >= 0x7f800000u)) + return as_special (x); // x=+inf, x=+nan + if (__glibc_likely (ax < 0x3c880000u)) + { // |x| < 0x1.1p-6 + if (__glibc_unlikely (ax < 0x33000000u)) + { // |x| < 0x1p-25 if (!ax) return x; return fmaf (x, -x, x); @@ -113,65 +110,43 @@ __log1pf (float x) 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])); + + z4 * (b[5] + z * (b[6] + z * b[7]))); double r = z + f; - if (__glibc_unlikely ((asuint64 (r) & 0xfffffffll) == 0)) + 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 - }; + uint64_t tp = asuint64 (z + 1.0); + int e = (tp >> 52) - 0x3ff; + uint64_t m52 = tp & (~0ull >> 12); + unsigned j = (tp >> (52 - 5)) & 31; + double xd = asdouble (m52 | UINT64_C (0x3ff) << 52); + z = xd * x0[j] - 1; // z is exact for x<0x1.0cp+30 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; + r = (ln2 * e + lix[j]) + + z * ((c[0] + z * c[1]) + z2 * (c[2] + z * c[3])); + const double eps = 2.1555e-11; + float ub = r, lb = r - eps; 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]))); + f = z2 + * ((b[1] + z * b[2]) + z2 * (b[3] + z * b[4]) + + z4 * (b[5] + z * (b[6] + z * b[7]))); + double lj = lix[j] - 0x1.7654p-37; // subtract the offset 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-6f) - return -0x1.271f0ep-6f - 0x1p-31f; - if (x == -0x1.3a415ep-5f) - return -0x1.407112p-5f + 0x1p-30f; - if (x == 0x1.fb035ap-2f) - return 0x1.9bddc2p-2f + 0x1p-27f; - tr += 64 * (rl + (Lh - tr)); - } - else if (rl + (Lh - tr) == 0.0) - { - if (x == 0x1.b7fd86p-4f) - return 0x1.a1ece2p-4f + 0x1p-29f; - if (x == -0x1.3a415ep-5f) - return -0x1.407112p-5f + 0x1p-30f; - if (x == 0x1.43c7e2p-6f) - return 0x1.409f80p-6f + 0x1p-31f; - } - ub = tr; + double Lh = ln2h * e, Ll = ln2l * e; + Ll += z; + double rh = Lh + lj, rl = ((Lh - rh) + lj) + (Ll + f); + float fh = rh + rl; + double Fl = (rh - (double) fh) + rl; + float fl = Fl, tfl = fl * 2.0f; + if ((fh + tfl) - fh == tfl) + fl += copysignf (0.5f, (float) (Fl - (double) fl)) * fabsf (fl); + ub = fh + fl; } return ub; }