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

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

• 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)
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
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]