For sufficiently small abs(x), gamma(x) is inaccurate.
Here is an example with x = 2^-256:
fpprec : 16 $ float(gamma(2b0^-256));
fpprec : 32 $ float(gamma(2b0^-256));
fpprec : 64 $ float(gamma(2b0^-256));
fpprec : 128 $ float(gamma(2b0^-256));
fpprec : 256 $ float(gamma(2b0^-256));
7.7675050341672e20
1.5574252028712945e37
-1.1009878111979649e69
1.157920892373162e77
1.157920892373162e77
Here is an example with x = -2^-256:
fpprec : 16 $ float(gamma(-2b0^-256));
fpprec : 32 $ float(gamma(-2b0^-256));
fpprec : 64 $ float(gamma(-2b0^-256));
fpprec : 128 $ float(gamma(-2b0^-256));
fpprec : 256 $ float(gamma(-2b0^-256));
7.7675050341672e20
1.5574252028712945e37
-1.1009878111979649e69
-1.157920892373162e77
-1.157920892373162e77
Even at 256 digits of precision, gamma(2b0^-256) has 72 incorrect trailing digits. It appears that, for any given precision, there is always a small enough input that will produce completely incorrect output. If we are using P decimal digits, the problem shows up at around abs(x) = 10^-P, at which point gamma(x) is very close to 1/x.
I am using Maxima 5.31.2 on a Windows machine.
Raymond Toy
2013-12-07
Thanks for the nice report. We can do a lot better than this by using the relationship
gamma(x) = %pi/sin(%pi*x)/gamma(1-x)
For such small x, gamma(1-x) is 1 and we get %pi/sin(%pi*x). For your example of 2^-256, we instantly get 1.157920892373162e+77 without any troubles.
Raymond Toy
2013-12-07
Raymond Toy
2013-12-07
Fixed in commit 0fb4df38c.
Dan Gildea
2013-12-15
The following test is failing with gcl and ecl:
closeto( gamma(2.0^-256 + 1d-200*%i), rectform(1/(2.0^-256 + 1d-200*%i)), 2d-15); true;
With cmucl:
Maxima branch_5_31_base_238_gf7e7ca9 http://maxima.sourceforge.net
using Lisp CMU Common Lisp 20c release-20c (20C Unicode)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) gamma(2.0^-256 + 1d-200*%i);
(%o1) 1.157920892373162e+77 - 1.34078079299426e-46 %i
(%i2) :lisp(trace gamma-lanczos)
(GAMMA-LANCZOS)
(%i2) gamma(2.0^-256 + 1d-200*%i);
0: (GAMMA-LANCZOS #C(8.636168555094445e-78 1.0e-200))
1: (GAMMA-LANCZOS #C(1.0 1.0e-200))
1: GAMMA-LANCZOS returned #C(1.0 -5.772156649015329e-201)
0: GAMMA-LANCZOS returned #C(1.157920892373162e77 -1.3407807929942597e-46)
(%o2) 1.157920892373162e+77 - 1.34078079299426e-46 %i
With ecl:
Maxima branch_5_31_base_238_gf7e7ca9 http://maxima.sourceforge.net
using Lisp ECL 13.5.1
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) :lisp(trace gamma-lanczos)
(GAMMA-LANCZOS)
(%i1) gamma(2.0^-256 + 1d-200*%i);
1> (GAMMA-LANCZOS #C(8.636168555094445e-78 1.e-200))
| 2> (GAMMA-LANCZOS #C(1.0 1.e-200))
| <2 (GAMMA-LANCZOS #C(0.9999999999999999 -5.772156649015328e-201))
<1 (GAMMA-LANCZOS #C(1.1579208923731618e77 -1.3407807929942595e-46))
(%o1) 1.1579208923731618e+77 - 13.407807929942595e-47 %i
(%i2)
Raymond Toy
2013-12-15
Looks like the threshold is a little bit too tight for ecl and gcl.
Dan Gildea
2013-12-16
Switched test to relerror in commit:
[fbcee9552427fa12f7223cde4a8069e270dc62cb]