Here is a patch for tammaq for negative argument pr68686.

I know about depending on ports from upstream but this was done recently and this (tgammaq) was left out.

This patch is basically a one-liner.

I have test cases but libquadmath doesn't have a testsuite.

One test just shows alternating signs for -0.5Q, -1.5Q, -2.5Q, etc.

Another test verifies the Gamma reflection formula:

 tgamma(1-x) * tgamma(x) - pi / sinq(pi*x) is tiny.

Builds on x86_64-linux and tests correctly offline.


OK?


2017-11-10  Edward Smith-Rowland  <3dw...@verizon.net>

        PR libquadmath/68686
        * math/tgammaq.c: Correct sign for negative argument.
diff --git a/libquadmath/math/tgammaq.c b/libquadmath/math/tgammaq.c
index a07d583..3080094 100644
--- a/libquadmath/math/tgammaq.c
+++ b/libquadmath/math/tgammaq.c
@@ -47,7 +47,9 @@ tgammaq (__float128 x)
     /* x == -Inf.  According to ISO this is NaN.  */
     return x - x;
 
-  /* XXX FIXME.  */
   res = expq (lgammaq (x));
-  return signbitq (x) ? -res : res;
+  if (x > 0.0Q || ((int)(-x) & 1) == 1)
+    return res;
+  else
+    return -res;
 }
#include <quadmath.h>
#include <assert.h>

void
test_alt_signs()
{
  __float128 result = tgammaq(-1.5Q);
  assert(result > 0.0Q);

  result = tgammaq(-2.5Q);
  assert(result < 0.0Q);

  result = tgammaq(-3.5Q);
  assert(result > 0.0Q);

  result = tgammaq(-4.5Q);
  assert(result < 0.0Q);
}

/*
 * Return |\Gamma(x) \Gamma(1 - x) - \frac{\pi}{sin(\pi x)}|
 */
__float128
abs_delta(__float128 x)
{
  return fabsq(tgammaq(x) * tgammaq(1.0Q - x) - M_PIq / sinq(M_PIq * x));
}

/*
 * Test reflection: \Gamma(x) \Gamma(1 - x) = \frac{\pi}{sin(\pi x)}
 */
void
test_reflection()
{
  assert(abs_delta(+1.5Q) < 100 * FLT128_EPSILON);
  assert(abs_delta(+0.5Q) < 100 * FLT128_EPSILON);
  assert(abs_delta(-0.5Q) < 100 * FLT128_EPSILON);
  assert(abs_delta(-1.5Q) < 100 * FLT128_EPSILON);
  assert(abs_delta(-2.5Q) < 100 * FLT128_EPSILON);
}

int
main()
{
  test_alt_signs();

  test_reflection();

  return 0;
}

Reply via email to