Improve the accuracy of tgamma for negative inputs (BZ 26983)
Commit Message
---
sysdeps/ieee754/dbl-64/e_gamma_r.c | 12 ++++++++++--
1 file changed, 10 insertions(+), 2 deletions(-)
Comments
On Fri, 26 Feb 2021, Paul Zimmermann wrote:
> ---
> sysdeps/ieee754/dbl-64/e_gamma_r.c | 12 ++++++++++--
> 1 file changed, 10 insertions(+), 2 deletions(-)
Please include a longer commit message explaining the change made and why
it's an appropriate fix for the issue seen.
> @@ -188,8 +189,15 @@ __ieee754_gamma_r (double x, int *signgamp)
> ? __sin (M_PI * frac)
> : __cos (M_PI * (0.5 - frac)));
> int exp2_adj;
> - double tret = M_PI / (-x * sinpix
> - * gamma_positive (-x, &exp2_adj));
> + double h1, l1, h2, l2;
> + mul_split (&h1, &l1, sinpix, gamma_positive (-x, &exp2_adj));
> + /* sinpix*gamma_positive(.) = h1 + l1 */
> + mul_split (&h2, &l2, h1, x);
> + /* h1*x = h2 + l2 */
> + /* (h1 + l1) * x = h1*x + l1*x = h2 + l2 + l1*x */
> + l2 += l1 * x;
> + h2 += l2;
> + double tret = M_PI / -h2;
Note that glibc uses tabs for the initial multiple-of-8-columns of
indentation; this is inserting code using spaces instead.
I'd expect similar changes for ldbl-96 and ldbl-128 at least, in the
absence of some reason it makes sense for the implementations to diverge.
The test input that justifies this change as being needed should also be
added to auto-libm-test-in by the patch (or XFAILs for that test input
removed if already present), with any corresponding libm-test-ulps updates
also being included for at least one architecture on which the patch was
tested. Likewise the test input you found giving too-large ulps for
binary128, unless the binary128 problem is still present after making a
corresponding fix to the implementation for ldbl-128.
On Feb 26 2021, Paul Zimmermann wrote:
> @@ -188,8 +189,15 @@ __ieee754_gamma_r (double x, int *signgamp)
> ? __sin (M_PI * frac)
> : __cos (M_PI * (0.5 - frac)));
> int exp2_adj;
> - double tret = M_PI / (-x * sinpix
> - * gamma_positive (-x, &exp2_adj));
> + double h1, l1, h2, l2;
> + mul_split (&h1, &l1, sinpix, gamma_positive (-x, &exp2_adj));
> + /* sinpix*gamma_positive(.) = h1 + l1 */
> + mul_split (&h2, &l2, h1, x);
> + /* h1*x = h2 + l2 */
> + /* (h1 + l1) * x = h1*x + l1*x = h2 + l2 + l1*x */
> + l2 += l1 * x;
> + h2 += l2;
> + double tret = M_PI / -h2;
There is a mix of tab and space indentation. Please always use tab
indentation.
Andreas.
@@ -24,6 +24,7 @@
#include <math-underflow.h>
#include <float.h>
#include <libm-alias-finite.h>
+#include <mul_split.h>
/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
approximation to gamma function. */
@@ -188,8 +189,15 @@ __ieee754_gamma_r (double x, int *signgamp)
? __sin (M_PI * frac)
: __cos (M_PI * (0.5 - frac)));
int exp2_adj;
- double tret = M_PI / (-x * sinpix
- * gamma_positive (-x, &exp2_adj));
+ double h1, l1, h2, l2;
+ mul_split (&h1, &l1, sinpix, gamma_positive (-x, &exp2_adj));
+ /* sinpix*gamma_positive(.) = h1 + l1 */
+ mul_split (&h2, &l2, h1, x);
+ /* h1*x = h2 + l2 */
+ /* (h1 + l1) * x = h1*x + l1*x = h2 + l2 + l1*x */
+ l2 += l1 * x;
+ h2 += l2;
+ double tret = M_PI / -h2;
ret = __scalbn (tret, -exp2_adj);
math_check_force_underflow_nonneg (ret);
}