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 <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
> #include
>
> 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
