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 <Paul.Zimmermann@inria.fr>

diff --git a/SHARED-FILES b/SHARED-FILES
index e942ae388ec02ad252c94c3ecc6acebbeda96e93..7550bc8b0b99d3245ba864e8fe1a5c64c4fee712 100644 (file)
--- 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 2001baa60558013988f8d1846990b1da11019aeb..4fd72b3ab6dc38cf42918e9b158572d3b8c35152 100644 (file)
--- 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 f7d3b35e6d4465c45ea6512fcd9c9cf6b695cdd3..1b5326a2d14e69ce378eb80fbe430b8027cf0022 100644 (file)
--- 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 f03c60df341600b280515068341bc94c332ce130..4cd5fa10aa13c0f454289a785afd0c59201f6d25 100644 (file)
--- 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;
     }