Fix tgamma missing underflows (bug 18951).

Similar to various other bugs in this area, tgamma functions can fail to raise the underflow exception when the result is tiny and inexact but one or more low bits of the intermediate result that is scaled down are zero. This patch forces the exception in a similar way to previous fixes. Tested for x86_64, x86, mips64 and powerpc. [BZ #18951] * sysdeps/ieee754/dbl-64/e_gamma_r.c (__ieee754_gamma_r): Force underflow exception for small results. * sysdeps/ieee754/flt-32/e_gammaf_r.c (__ieee754_gammaf_r): Likewise. * sysdeps/ieee754/ldbl-128/e_gammal_r.c (__ieee754_gammal_r): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c (__ieee754_gammal_r): Likewise. * sysdeps/ieee754/ldbl-96/e_gammal_r.c (__ieee754_gammal_r): Likewise. * math/auto-libm-test-in: Add more tests of tgamma. * math/auto-libm-test-out: Regenerated.
git-mirror · Sep 17, 2015 · 46f74e1 · 46f74e1
1 parent e67dc1b
commit 46f74e1
Show file tree

Hide file tree

Showing 9 changed files with 1,084 additions and 2 deletions.
diff --git a/ChangeLog b/ChangeLog
@@ -1,3 +1,19 @@
+2015-09-17  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #18951]
+	* sysdeps/ieee754/dbl-64/e_gamma_r.c (__ieee754_gamma_r): Force
+	underflow exception for small results.
+	* sysdeps/ieee754/flt-32/e_gammaf_r.c (__ieee754_gammaf_r):
+	Likewise.
+	* sysdeps/ieee754/ldbl-128/e_gammal_r.c (__ieee754_gammal_r):
+	Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c (__ieee754_gammal_r):
+	Likewise.
+	* sysdeps/ieee754/ldbl-96/e_gammal_r.c (__ieee754_gammal_r):
+	Likewise.
+	* math/auto-libm-test-in: Add more tests of tgamma.
+	* math/auto-libm-test-out: Regenerated.
+
 2015-09-17  Andreas Schwab  <schwab@suse.de>
 
 	* sysdeps/unix/sysv/linux/sparc/sparc64/semctl.c: Readd.

diff --git a/NEWS b/NEWS
@@ -14,8 +14,8 @@ Version 2.23
   17244, 17787, 17905, 18084, 18086, 18240, 18265, 18370, 18421, 18480,
   18525, 18595, 18610, 18618, 18647, 18661, 18674, 18675, 18681, 18757,
   18778, 18781, 18787, 18789, 18790, 18795, 18796, 18820, 18823, 18824,
-  18857, 18863, 18870, 18872, 18873, 18875, 18887, 18921, 18952, 18961,
-  18966, 18967, 18977.
+  18857, 18863, 18870, 18872, 18873, 18875, 18887, 18921, 18951, 18952,
+  18961, 18966, 18967, 18977.
 
 * The obsolete header <regexp.h> has been removed.  Programs that require
   this header must be updated to use <regex.h> instead.

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
@@ -3855,6 +3855,47 @@ tgamma 1e3
 tgamma -100000.5
 tgamma max
 
+tgamma -0x22.30p0
+tgamma -0x22.31p0
+tgamma -0x22.32p0
+tgamma -0x22.33p0
+tgamma -0x22.34p0
+tgamma -0x22.35p0
+tgamma -0x22.36p0
+tgamma -0x22.37p0
+tgamma -0xa3.70p0
+tgamma -0xa3.71p0
+tgamma -0xa3.72p0
+tgamma -0xa3.73p0
+tgamma -0xa3.74p0
+tgamma -0xa3.75p0
+tgamma -0xa3.76p0
+tgamma -0xa3.77p0
+tgamma -0xab.0d0p0
+tgamma -0xab.0d1p0
+tgamma -0xab.0d2p0
+tgamma -0xab.0d3p0
+tgamma -0xab.0d4p0
+tgamma -0xab.0d5p0
+tgamma -0xab.0d6p0
+tgamma -0xab.0d7p0
+tgamma -0x6db.030p0
+tgamma -0x6db.031p0
+tgamma -0x6db.032p0
+tgamma -0x6db.033p0
+tgamma -0x6db.034p0
+tgamma -0x6db.035p0
+tgamma -0x6db.036p0
+tgamma -0x6db.037p0
+tgamma -0x6db.050p0
+tgamma -0x6db.051p0
+tgamma -0x6db.052p0
+tgamma -0x6db.053p0
+tgamma -0x6db.054p0
+tgamma -0x6db.055p0
+tgamma -0x6db.056p0
+tgamma -0x6db.057p0
+
 tgamma -0x3.06644cp+0
 tgamma -0x6.fe4636e0c5064p+0
 tgamma -0x7.a13d7a2945cd5718p+0

diff --git a/math/auto-libm-test-out b/math/auto-libm-test-out
diff --git a/sysdeps/ieee754/dbl-64/e_gamma_r.c b/sysdeps/ieee754/dbl-64/e_gamma_r.c
@@ -194,6 +194,11 @@ __ieee754_gamma_r (double x, int *signgamp)
 	      double tret = M_PI / (-x * sinpix
 				    * gamma_positive (-x, &exp2_adj));
 	      ret = __scalbn (tret, -exp2_adj);
+	      if (ret < DBL_MIN)
+		{
+		  double force_underflow = ret * ret;
+		  math_force_eval (force_underflow);
+		}
 	    }
 	}
     }

diff --git a/sysdeps/ieee754/flt-32/e_gammaf_r.c b/sysdeps/ieee754/flt-32/e_gammaf_r.c
@@ -186,6 +186,11 @@ __ieee754_gammaf_r (float x, int *signgamp)
 	      float tret = (float) M_PI / (-x * sinpix
 					   * gammaf_positive (-x, &exp2_adj));
 	      ret = __scalbnf (tret, -exp2_adj);
+	      if (ret < FLT_MIN)
+		{
+		  float force_underflow = ret * ret;
+		  math_force_eval (force_underflow);
+		}
 	    }
 	}
     }

diff --git a/sysdeps/ieee754/ldbl-128/e_gammal_r.c b/sysdeps/ieee754/ldbl-128/e_gammal_r.c
@@ -194,6 +194,11 @@ __ieee754_gammal_r (long double x, int *signgamp)
 	      ret = M_PIl / (-x * sinpix
 			     * gammal_positive (-x, &exp2_adj));
 	      ret = __scalbnl (ret, -exp2_adj);
+	      if (ret < LDBL_MIN)
+		{
+		  long double force_underflow = ret * ret;
+		  math_force_eval (force_underflow);
+		}
 	    }
 	}
     }

diff --git a/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c b/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c
@@ -194,6 +194,11 @@ __ieee754_gammal_r (long double x, int *signgamp)
 	      ret = M_PIl / (-x * sinpix
 			     * gammal_positive (-x, &exp2_adj));
 	      ret = __scalbnl (ret, -exp2_adj);
+	      if (ret < LDBL_MIN)
+		{
+		  long double force_underflow = ret * ret;
+		  math_force_eval (force_underflow);
+		}
 	    }
 	}
     }

diff --git a/sysdeps/ieee754/ldbl-96/e_gammal_r.c b/sysdeps/ieee754/ldbl-96/e_gammal_r.c
@@ -186,6 +186,11 @@ __ieee754_gammal_r (long double x, int *signgamp)
 	      ret = M_PIl / (-x * sinpix
 			     * gammal_positive (-x, &exp2_adj));
 	      ret = __scalbnl (ret, -exp2_adj);
+	      if (ret < LDBL_MIN)
+		{
+		  long double force_underflow = ret * ret;
+		  math_force_eval (force_underflow);
+		}
 	    }
 	}
     }