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.