[COMMITTED] math: Sync tgamma with CORE-MATH

Adhemerval Zanella adhemerval.zanella@linaro.org
Mon Sep 21 12:39:08 GMT 2026


Import the following changes from CORE-MATH (up to revision c537979a):

  - The error bounds of two fast paths were too tight to cover all the
    hard-to-round cases found by exhaustive search, so the rounding test
    wrongly accepted the fast-path result.

  - The conversion of floor(x) to int64_t is undefined for x below
    INT64_MIN, so k is now computed only on the paths that use it, and
    clamped on the negative non-integer path.

  - The comparisons against 0.0f on the integer path are made against 0.0
    (x is a double).

Checked on aarch64-linux-gnu.
---
 SHARED-FILES                       |  2 +-
 sysdeps/ieee754/dbl-64/e_gamma_r.c | 14 ++++++++------
 2 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/SHARED-FILES b/SHARED-FILES
index b856ef9a87..86ad228372 100644
--- a/SHARED-FILES
+++ b/SHARED-FILES
@@ -269,7 +269,7 @@ core-math:
   sysdeps/ieee754/dbl-64/e_atanh.c
   # src/binary64/cosh/cosh.c, revision 5d0c89d5
   sysdeps/ieee754/dbl-64/e_cosh.c
-  # src/binary64/tgamma/tgamma.c, revision 0f185e23
+  # src/binary64/tgamma/tgamma.c, revision c537979a
   sysdeps/ieee754/dbl-64/e_gamma_r.c
   # src/binary64/lgamma/lgamma.c, revision c537979a
   sysdeps/ieee754/dbl-64/e_lgamma_r.c
diff --git a/sysdeps/ieee754/dbl-64/e_gamma_r.c b/sysdeps/ieee754/dbl-64/e_gamma_r.c
index 997c63a7df..a94381cede 100644
--- a/sysdeps/ieee754/dbl-64/e_gamma_r.c
+++ b/sysdeps/ieee754/dbl-64/e_gamma_r.c
@@ -3,7 +3,7 @@
 Copyright (c) 2024-2025 Alexei Sibidanov <sibid@uvic.ca>.
 
 The original version of this file was copied from the CORE-MATH
-project (file src/binary64/tgamma/tgamma.c, revision 0f185e23).
+project (file src/binary64/tgamma/tgamma.c, revision c537979a).
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
@@ -772,13 +772,14 @@ __ieee754_gamma_r (double x, int *signgamp)
   double fx = floor (x);
   /* compute k only after the overflow check, otherwise the cast to integer
      might overflow */
-  int64_t k = fx;
+  int64_t k;
   if (__glibc_unlikely (fx == x))
     { /* x is integer */
-      if (x == 0.0f)
+      if (x == 0.0)
 	return __math_divzero (0);
-      if (x < 0.0f)
+      if (x < 0.0)
 	return __math_invalid (0);
+      k = fx;
       double t0h = 1, t0l = 0, x0 = 1;
       for (int i = 1; i < k; i++, x0 += 1.0)
 	t0h = mulddd2 (x0, t0h, t0l, &t0l);
@@ -789,6 +790,7 @@ __ieee754_gamma_r (double x, int *signgamp)
     { /* negative non-integer */
       /* For x <= -184, x non-integer, |gamma(x)| < 2^-1078.  */
       static const double sgn[2] = { 0x1p-1022, -0x1p-1022 };
+      k = fx < (double) INT64_MIN ? INT64_MIN : fx;
       return __math_erange (0x1p-1022 * sgn[k & 1]);
     }
 
@@ -810,7 +812,7 @@ __ieee754_gamma_r (double x, int *signgamp)
 	  rh = -rh;
 	  rl = -rl;
 	}
-      double eps = rh * (0x1.eb2049057bc61p-68 - x * 0x1.61019f74442b7p-73);
+      double eps = rh * (0x1.48ad3caefe28p-67 - x * 0x1.61019f74442b7p-73);
       uint64_t th;
       if (__glibc_likely (ip >= -170))
 	{
@@ -855,7 +857,7 @@ __ieee754_gamma_r (double x, int *signgamp)
       double ll = 0, lh = as_lgamma_asym (x, &ll);
       int e;
       lh = as_expd (lh, &ll, &e);
-      double eps = lh * (0x1.2e3b40a0e9b4fp-69 + x * 0x1.6aad80c11872cp-73);
+      double eps = lh * (0x1.2e3b40a0e9b4fp-69 + x * 0x1.bce27bce24d11p-73);
       double ub = lh + (ll + eps), lb = lh + (ll - eps);
       if (ub != lb)
 	return as_tgamma_accurate (x);
-- 
2.47.3



More information about the Libc-alpha mailing list