[COMMITTED] math: Sync lgammaf with CORE-MATH

[Paul Zimmermann]

       Hi,

while re-running exhaustive checks with CORE-MATH test suite, I hit:

zimmerma@thym:~/svn/core-math$ CORE_MATH_CHECK_STD=true CORE_MATH_LAUNCHER="/tmp/lib/ld-linux-x86-64.so.2 --library-path /tmp/lib:/usr/lib/x86_64-linux-gnu" LDFLAGS="-L /tmp/lib" ./check.sh lgammaf
Running exhaustive check in --rndn mode...
Error, signgam wrong for x=-0x1.ecp+6
expected 1, got -1

Is signgam defined for negative integers? The reference manual is not very
clear about this:

       The  sign  of  the
       Gamma  function  is returned in the external integer signgam declared in
       <math.h>.  It is 1 when the Gamma function is positive or zero, -1  when
       it is negative.

I see on https://pubs.opengroup.org/onlinepubs/9799919799/functions/lgamma.html:

If \(x\) is NaN, -Inf, or a negative integer, the value of signgam is unspecified.

Should we add this in the manual?

Paul

> From: Adhemerval Zanella <adhemerval.zanella@linaro.org>
> Date: Thu, 19 Mar 2026 13:53:20 -0300
> 
> 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
> 
>