On 11/3/18 10:09 PM, Jeff Law wrote:
On 10/23/18 7:45 PM, Ed Smith-Rowland wrote:
Greetings,
This is an almost trivial patch to get the correct sign for tgammaq.
I don't have a testcase as I don't know where to put one.
OK?
Ed Smith-Rowland
tgammaq.CL
2018-10-24 Edward Smith-Rowland <3dw...@verizon.net>
PR libquadmath/68686
* math/tgammaq.c: Correct sign for negative argument.
I don't have the relevant background to evaluate this for correctness.
Can you refer back to any kind of documentation which indicates what the
sign of the return value ought to be?
Alternately, if you can point to the relevant code in glibc that handles
the resultant sign, that'd be useful too.
Note that Joseph's follow-up doesn't touch on the gamma problem AFAICT,
but instead touches on the larger issues around trying to keep the
quadmath implementations between glibc and gcc more in sync.
Jeff
Thank you for (re)considering this.The maths here can be read from a
very good NIST website:
For the functional formulas referred to below - https://dlmf.nist.gov/5.5
For cool pictures - https://dlmf.nist.gov/5.3#i
(Aside: the reciprocal gamma function is entire - possessing no poles or
other analytic complications anywhere in the complex plane. A rgamma
function might be a nice addition to the C family and to TS 18661-1 for
that matter.)
For a table of extrema (Table 5.4.1) - https://dlmf.nist.gov/5.4#iii
These would be nice tests for accuracy as well as sign.
TL;DR:
======
Either follow the factorial-like recursion Gamma(x+1) = x Gamma(x) [DLMF
5.5.1] backwards:
Gamma(x-1) = Gamma(x) / (x-1)
Given that Gamma(x) is positive for x > 0 this will give alternating
negative signs if you start with, say, x=1/2 and keep going.
A cooler formula is [DLMF 5.5.3]:
Gamma(x) Gamma(1-x) = pi / sin(pi x)
Start with x = 3/2 and continue with higher odd half integers to see the
sign alternation.
GLibC
=====
I looked in glibc. Unfortunately, I see how they have the same mistake:
glibc/math/w_tgammal_compat.c:
long double
__tgammal(long double x)
{
int local_signgam;
long double y = __ieee754_gammal_r(x,&local_signgam);
...
return local_signgam < 0 ? - y : y;
}
I'm very sure this is where tgammaq came from.
Ditto for glibc/math/w_tgamma_compat.c and glibc/math/w_tgammaf_compat.c.
This fix will need to be done upstream.
Ed
