Ed, gfortran uses libquadmath for support of its REAL(16) intrinsic subprograms. I checked the list of intrinsics and tgamma is current not in the list. I don't have time to try Fortran's C interop feature to see if an ordinary user can access __float128. While your patch is probably useful for GCC, I doubr that gfortran in its current state will use tgammal. You'll need either Jakub or Joseph (jsm28) to give an OK.
-- steve On Mon, Nov 13, 2017 at 01:27:42PM -0500, Ed Smith-Rowland wrote: > 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 <[email protected]> > > 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; > } -- Steve 20170425 https://www.youtube.com/watch?v=VWUpyCsUKR4 20161221 https://www.youtube.com/watch?v=IbCHE-hONow
