]> git.ipfire.org Git - thirdparty/glibc.git/commitdiff
math: Fix non-portability in the computation of signgam in lgammaf
authorVincent Lefevre <vincent@vinc17.net>
Fri, 22 Nov 2024 16:54:53 +0000 (13:54 -0300)
committerAdhemerval Zanella <adhemerval.zanella@linaro.org>
Mon, 25 Nov 2024 12:20:47 +0000 (09:20 -0300)
The k>>31 in signgam = 1 - (((k&(k>>31))&1)<<1); is not portable:

* The ISO C standard says "If E1 has a signed type and a negative
  value, the resulting value is implementation-defined." (this is
  still in C23).
* If the int type is larger than 32 bits (e.g. a 64-bit type),
  then k = INT_MAX; line 144 will make k>>31 put 1 in bit 0
  (thus signgam will be -1) while 0 is expected.

Moreover, instead of the fx >= 0x1p31f condition, testing fx >= 0
is probably better for 2 reasons:

The signgam expression has more or less a condition on the sign
of fx (the goal of k>>31, which can be dropped with this new
condition). Since fx ≥ 0 should be the most common case, one can
get signgam directly in this case (value 1). And this simplifies
the expression for the other case (fx < 0).

This new condition may be easier/faster to test on the processor
(e.g. by avoiding a load of a constant from the memory).

This is commit d41459c731865516318f813cf4c966dafa0eecbf from CORE-MATH.

Checked on x86_64-linux-gnu.

sysdeps/ieee754/flt-32/e_lgammaf_r.c

index cb6551305605ed87c9e038068ff437461e150783..447376dc552039dc36afa4509e26b159d09cf51d 100644 (file)
@@ -181,12 +181,11 @@ __ieee754_lgammaf_r (float x, int *signgamp)
      Note that for a binary32 |x| >= 2^23, x is necessarily an integer,
      and we already dealed with negative integers, thus now:
      -2^23 < x < +Inf and x is not a negative integer nor 0, 1, 2. */
-  int k;
-  if (__builtin_expect (fx >= 0x1p31f, 0))
-    k = INT_MAX;
+  if (__glibc_unlikely (fx >= 0))
+    *signgamp = 1;
   else
-    k = fx;
-  *signgamp = 1 - (((k & (k >> 31)) & 1) << 1);
+    /* gamma(x) is negative in (-2n-1,-2n), thus when fx is odd.  */
+    *signgamp = 1 - ((((int) fx) & 1) << 1);
 
   double z = ax, f;
   if (__glibc_unlikely (ax < 0x1.52p-1f))