[COMMITTED] math: Sync lgammaf with CORE-MATH

Adhemerval Zanella adhemerval.zanella@linaro.org
Thu Mar 19 16:53:20 GMT 2026


It removes some unnecessary corner-case checks and uses a slightly
different binary algorithm for the hard-case database binary search.

Checked on aarch64-linux-gnu, arm-linux-gnueabihf,
powerpc64le-linux-gnu, i686-linux-gnu, and x86_64-linux-gnu.
---
 SHARED-FILES                         |  2 +-
 sysdeps/ieee754/flt-32/e_lgammaf_r.c | 36 ++++++++++------------------
 2 files changed, 14 insertions(+), 24 deletions(-)

diff --git a/SHARED-FILES b/SHARED-FILES
index c3e8e0f5ae..5a676e1a48 100644
--- a/SHARED-FILES
+++ b/SHARED-FILES
@@ -272,7 +272,7 @@ core-math:
   sysdeps/ieee754/flt-32/e_coshf.c
   # src/binary32/tgamma/tgammaf.c, revision 8ea8ea35
   sysdeps/ieee754/flt-32/e_gammaf_r.c
-  # src/binary32/lgamma/lgammaf.c, revision bc385c2
+  # src/binary32/lgamma/lgammaf.c, revision 8ea8ea35
   sysdeps/ieee754/flt-32/e_lgammaf_r.c
   # src/binary32/log10/log10f.c, revision ebff4c43
   sysdeps/ieee754/flt-32/e_log10f.c
diff --git a/sysdeps/ieee754/flt-32/e_lgammaf_r.c b/sysdeps/ieee754/flt-32/e_lgammaf_r.c
index 2fb4784b59..7b947cd979 100644
--- a/sysdeps/ieee754/flt-32/e_lgammaf_r.c
+++ b/sysdeps/ieee754/flt-32/e_lgammaf_r.c
@@ -1,10 +1,10 @@
 /* Correctly-rounded logarithm of the absolute value of the gamma function
    for binary32 value.
 
-Copyright (c) 2023, 2024 Alexei Sibidanov.
+Copyright (c) 2023-2026 Alexei Sibidanov.
 
 This file is part of the CORE-MATH project
-project (file src/binary32/lgamma/lgammaf.c, revision bc385c2).
+project (file src/binary32/lgamma/lgammaf.c, revision 8ea8ea35).
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
@@ -159,7 +159,7 @@ __lgammaf_r (float x, int *signgamp)
       return x + x; /* nan */
     }
   if (__glibc_unlikely (fx == x))
-    {
+    { /* x integer */
       if (x <= 0.0f)
 	{
 	  *signgamp = asuint (x) >> 31 ? -1 : 1;
@@ -204,15 +204,11 @@ __lgammaf_r (float x, int *signgamp)
       f = (c0 * s) * as_r8 (s, rn) / as_r8 (s, rd) - as_ln (z);
     }
   else
-    {
+    { /* |x| >= 0x1.52p-1 */
       if (ax > 0x1.afc1ap+1f)
 	{
-	  if (__glibc_unlikely (x > 0x1.895f1cp+121f))
+	  if (__glibc_unlikely (x >= 0x1.895f1cp+121f))
 	    return __math_oflowf (0);
-
-	  /* |x|>=2**23, must be -integer */
-	  if (__glibc_unlikely (x < 0.0f && ax > 0x1p+23f))
-	    return __math_divzerof (0);
 	  double lz = as_ln (z);
 	  f = (z - 0.5) * (lz - 1) + 0x1.acfe390c97d69p-2;
 	  if (ax < 0x1.0p+20f)
@@ -270,12 +266,6 @@ __lgammaf_r (float x, int *signgamp)
 	      -0x1.3a6c8295b4445p-1, -0x1.da44e8b810024p-3,
 	      -0x1.9061e81c77e4ap-5
 	    };
-	  if (x < 0.0f)
-	    {
-	      int ni = floorf (-2 * x);
-	      if ((ni & 1) == 0 && ni == -2 * x)
-		return __math_divzerof (0);
-	    }
 	  const double c0 = 0x1.3cc0e6a0106b3p+2;
 	  static const double rd[] =
 	    {
@@ -351,15 +341,15 @@ __lgammaf_r (float x, int *signgamp)
   if (__glibc_unlikely (tl <= 31u))
     {
       t = asuint (x);
-      int a = 0, b = array_length (tb) - 1;
-      while (a < b)
-	{ /* Binary search.  */
-	  int m = (a + b) >> 1;
-	  uint32_t tbi = asuint (tb[m].x);
-	  if (t > tbi)
-	    a = m + 1;
+      int a = 0, b = array_length (tb);
+      /* invariant: t.u < tb[0].x.u or tb[a].x.u <= t.u < tb[b].x.u */
+      while (a + 1 < b)
+	{
+	  int i = (a + b) / 2;
+	  if (t < asuint (tb[i].x))
+	    b = i;
 	  else
-	    b = m;
+	    a = i;
 	}
       if (t == asuint (tb[a].x))
 	return tb[a].f + tb[a].df;
-- 
2.43.0



More information about the Libc-alpha mailing list