Re: [PATCH libquadmath, pr68686] tgammaq(x) is always negative for noninteger x < 0

[Steve Kargl <sgk at troutmask dot apl dot washington dot edu>]

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  <3dw4rd@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;
> }


-- 
Steve
20170425 https://www.youtube.com/watch?v=VWUpyCsUKR4
20161221 https://www.youtube.com/watch?v=IbCHE-hONow