There are effectively two versions of the tgamma
function internally: a fully generic version that is slow, but reasonably
accurate, and a much more efficient approximation that is used where the
number of digits in the significand of T correspond to a certain Lanczos
approximation. In practice any built in floating point type you will
encounter has an appropriate Lanczos
approximation defined for it. It is also possible, given enough machine
time, to generate further Lanczos approximation's
using the program libs/math/tools/lanczos_generator.cpp.

The return type of this function is computed using the result
type calculation rules: the result is double
when T is an integer type, and T otherwise.

Returns tgamma(dz+1)-1.
Internally the implementation does not make use of the addition and subtraction
implied by the definition, leading to accurate results even for very small
dz. However, the implementation
is capped to either 35 digit accuracy, or to the precision of the Lanczos
approximation associated with type T, whichever is more accurate.

The return type of this function is computed using the result
type calculation rules: the result is double
when T is an integer type, and T otherwise.

The following table shows the peak errors (in units of epsilon) found on
various platforms with various floating point types, along with comparisons
to the GSL-1.9, GNU C Lib, HP-UX
C Library and Cephes
libraries. Unless otherwise specified any floating point type that is narrower
than the one shown will have effectively
zero error.

The gamma is relatively easy to test: factorials and half-integer factorials
can be calculated exactly by other means and compared with the gamma function.
In addition, some accuracy tests in known tricky areas were computed at high
precision using the generic version of this function.

The function tgamma1pm1 is
tested against values calculated very naively using the formula tgamma(1+dz)-1 with a
lanczos approximation accurate to around 100 decimal digits.

The generic version of the tgamma
function is implemented Sterling's approximation for lgamma for large z:

Following exponentiation, downward recursion is then used for small values
of z.

For types of known precision the Lanczos
approximation is used, a traits class boost::math::lanczos::lanczos_traits
maps type T to an appropriate approximation.

For z in the range -20 < z < 1 then recursion is used to shift to z
> 1 via:

For very small z, this helps to preserve the identity:

For z < -20 the reflection formula:

is used. Particular care has to be taken to evaluate the z * sin(π *
z) part: a special routine is used to reduce z prior to multiplying
by π to ensure that the result in is the range [0, π/2]. Without this an excessive
amount of error occurs in this region (which is hard enough already, as the
rate of change near a negative pole is exceptionally
high).

Finally if the argument is a small integer then table lookup of the factorial
is used.

The function tgamma1pm1 is
implemented using rational approximations devised
by JM in the region -0.5<dz<2.
These are the same approximations (and internal routines) that are used for
lgamma, and so aren't
detailed further here. The result of the approximation is log(tgamma(dz+1)) which can
fed into expm1 to give the
desired result. Outside the range -0.5<dz<2
then the naive formula tgamma1pm1(dz)=tgamma(dz+1)-1
can be used directly.