#2668 Bigfloat Gamma Inaccurate for Small Inputs

None
closed
Raymond Toy
5
2013-12-16
2013-12-06
Dane Vandeputte
No

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.

Discussion

  • Raymond Toy
    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
    Raymond Toy
    2013-12-07

    • labels: --> Lisp Core - Floating point
    • status: open --> closed
    • assigned_to: Raymond Toy
     
  • Raymond Toy
    Raymond Toy
    2013-12-07

    Fixed in commit 0fb4df38c.

     
  • Dan Gildea
    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
    Raymond Toy
    2013-12-15

    Looks like the threshold is a little bit too tight for ecl and gcl.