Fix tgamma missing underflows (bug 18951) [committed]
Commit Message
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. Committed.
(auto-libm-test-out diffs omitted below.)
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.
@@ -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
@@ -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);
+ }
}
}
}
@@ -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);
+ }
}
}
}
@@ -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);
+ }
}
}
}
@@ -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);
+ }
}
}
}
@@ -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);
+ }
}
}
}