[COMMITTED] math: Sync lgammaf with CORE-MATH

[Adhemerval Zanella Netto]

On 20/03/26 07:42, Paul Zimmermann wrote:
>        Hi Adhemerval,
> 
> with the CORE-MATH test suite, I get:
> 
> 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...
> all ok
> Running exhaustive check in --rndz mode...
> Spurious overflow exception for x=0x1.895f1cp+121 (y=0x1.fffffep+127)
> 
> Please can you check?

It seems to be a pre-existing issue, on 5c28a3cf9f (pre-sync) if I add
0x1.895f1cp+121  on math/auto-libm-test-in and regenerating 
math/auto-libm-test-out-lgamma (as below) I see:

$ math/test-float-lgamma
testing float (without inline functions)
Failure: lgamma (0x3.12be38p+120): errno set to 0, expected 34 (ERANGE)
Failure: lgamma_upward (0x3.12be38p+120): errno set to 0, expected 34 (ERANGE)

Test suite completed:
  Maximum error found of `0' ulp
  8 max error test cases,
  7368 input tests,
  - with 29488 tests for exception flags,
  - with 3688 tests for errno executed.
  2 errors occurred.

I forgot to sync the some ERANGE fixes like 20d51f2ee1e945cdc910aa3fbd40c7d5959ff3d0.
I will add the following patch:

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 4c3cd23ffd..2574957ada 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -6994,6 +6994,7 @@ lgamma -0x3.ec4298p+0
 lgamma 0x1.ecf3fep-73
 lgamma 0x1.58ace8p+112
 lgamma -0x1.efc2a2p+14
+lgamma 0x1.895f1cp+121

 # Values +/- 10ulp from overflow threshold.  (Values very close to
 # overflow threshold produce results very close of that threshold,
diff --git a/math/auto-libm-test-out-lgamma b/math/auto-libm-test-out-lgamma
index d27c186639..fed3f86c6c 100644
--- a/math/auto-libm-test-out-lgamma
+++ b/math/auto-libm-test-out-lgamma
@@ -2301,6 +2301,31 @@ lgamma -0x1.efc2a2p+14
 = lgamma tonearest ibm128 -0x7.bf0a88p+12 : -0x4.88b6f00000005971e29c3b8dd8p+16 -1 : inexact-ok
 = lgamma towardzero ibm128 -0x7.bf0a88p+12 : -0x4.88b6f00000005971e29c3b8dd8p+16 -1 : inexact-ok
 = lgamma upward ibm128 -0x7.bf0a88p+12 : -0x4.88b6f00000005971e29c3b8dd8p+16 -1 : inexact-ok
+lgamma 0x1.895f1cp+121
+= lgamma downward binary32 0x3.12be38p+120 : 0xf.fffffp+124 1 : inexact-ok
+= lgamma tonearest binary32 0x3.12be38p+120 : plus_infty 1 : inexact-ok overflow errno-erange
+= lgamma towardzero binary32 0x3.12be38p+120 : 0xf.fffffp+124 1 : inexact-ok
+= lgamma upward binary32 0x3.12be38p+120 : plus_infty 1 : inexact-ok overflow errno-erange
+= lgamma downward binary64 0x3.12be38p+120 : 0xf.fffffbf5af3fp+124 1 : inexact-ok
+= lgamma tonearest binary64 0x3.12be38p+120 : 0xf.fffffbf5af3f8p+124 1 : inexact-ok
+= lgamma towardzero binary64 0x3.12be38p+120 : 0xf.fffffbf5af3fp+124 1 : inexact-ok
+= lgamma upward binary64 0x3.12be38p+120 : 0xf.fffffbf5af3f8p+124 1 : inexact-ok
+= lgamma downward intel96 0x3.12be38p+120 : 0xf.fffffbf5af3f533p+124 1 : inexact-ok
+= lgamma tonearest intel96 0x3.12be38p+120 : 0xf.fffffbf5af3f534p+124 1 : inexact-ok
+= lgamma towardzero intel96 0x3.12be38p+120 : 0xf.fffffbf5af3f533p+124 1 : inexact-ok
+= lgamma upward intel96 0x3.12be38p+120 : 0xf.fffffbf5af3f534p+124 1 : inexact-ok
+= lgamma downward m68k96 0x3.12be38p+120 : 0xf.fffffbf5af3f533p+124 1 : inexact-ok
+= lgamma tonearest m68k96 0x3.12be38p+120 : 0xf.fffffbf5af3f534p+124 1 : inexact-ok
+= lgamma towardzero m68k96 0x3.12be38p+120 : 0xf.fffffbf5af3f533p+124 1 : inexact-ok
+= lgamma upward m68k96 0x3.12be38p+120 : 0xf.fffffbf5af3f534p+124 1 : inexact-ok
+= lgamma downward binary128 0x3.12be38p+120 : 0xf.fffffbf5af3f533cd835e6dc7ed8p+124 1 : inexact-ok
+= lgamma tonearest binary128 0x3.12be38p+120 : 0xf.fffffbf5af3f533cd835e6dc7ed8p+124 1 : inexact-ok
+= lgamma towardzero binary128 0x3.12be38p+120 : 0xf.fffffbf5af3f533cd835e6dc7ed8p+124 1 : inexact-ok
+= lgamma upward binary128 0x3.12be38p+120 : 0xf.fffffbf5af3f533cd835e6dc7eep+124 1 : inexact-ok
+= lgamma downward ibm128 0x3.12be38p+120 : 0xf.fffffbf5af3f533cd835e6dc7cp+124 1 : inexact-ok
+= lgamma tonearest ibm128 0x3.12be38p+120 : 0xf.fffffbf5af3f533cd835e6dc8p+124 1 : inexact-ok
+= lgamma towardzero ibm128 0x3.12be38p+120 : 0xf.fffffbf5af3f533cd835e6dc7cp+124 1 : inexact-ok
+= lgamma upward ibm128 0x3.12be38p+120 : 0xf.fffffbf5af3f533cd835e6dc8p+124 1 : inexact-ok
 lgamma 0x3.12be0cp+120
 = lgamma downward binary32 0x3.12be0cp+120 : 0xf.ffff1p+124 1 : inexact-ok
 = lgamma tonearest binary32 0x3.12be0cp+120 : 0xf.ffff1p+124 1 : inexact-ok
diff --git a/sysdeps/ieee754/flt-32/e_lgammaf_r.c b/sysdeps/ieee754/flt-32/e_lgammaf_r.c
index 7b947cd979..8efec19b48 100644
--- a/sysdeps/ieee754/flt-32/e_lgammaf_r.c
+++ b/sysdeps/ieee754/flt-32/e_lgammaf_r.c
@@ -208,7 +208,19 @@ __lgammaf_r (float x, int *signgamp)
       if (ax > 0x1.afc1ap+1f)
        {
          if (__glibc_unlikely (x >= 0x1.895f1cp+121f))
-           return __math_oflowf (0);
+           {
+             /* for x=0x1.895f1cp+121, lgamma(x) < 2^128, thus there is no
+                overflow for rounding towards zero or downwards.
+                The following expression overflows for x > 0x1.895f1cp+121
+                or x = 0x1.895f1cp+121 and rounding to nearest or away,
+                and does not overflow for x = 0x1.895f1cp+121 and rounding
+                towards zero */
+             float r = fmaf (x, 0x1.4d3398p+6f, 0x1.10f35ep+103f);
+             if (x > 0x1.895f1cp+121f || (x == 0x1.895f1cp+121f &&
+                 x * 5.0f >= 0x1.ebb6e4p+123))
+               return __math_oflowf (0);
+             return r;
+           }
          double lz = as_ln (z);
          f = (z - 0.5) * (lz - 1) + 0x1.acfe390c97d69p-2;
          if (ax < 0x1.0p+20f)