...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
#include <boost/math/special_functions/gamma.hpp>
namespace boost{ namespace math{ template <class T> calculatedresulttype tgamma(T z); template <class T, class Policy> calculatedresulttype tgamma(T z, const Policy&); template <class T> calculatedresulttype tgamma1pm1(T dz); template <class T, class Policy> calculatedresulttype tgamma1pm1(T dz, const Policy&); }} // namespaces
template <class T> calculatedresulttype tgamma(T z); template <class T, class Policy> calculatedresulttype tgamma(T z, const Policy&);
Returns the "true gamma" (hence name tgamma) of value z:
The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.
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.
template <class T> calculatedresulttype tgamma1pm1(T dz); template <class T, class Policy> calculatedresulttype tgamma1pm1(T dz, const Policy&);
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
.
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 final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.
The following table shows the peak errors (in units of epsilon) found on various platforms with various floating point types, along with comparisons to other common libraries. Unless otherwise specified any floating point type that is narrower than the one shown will have effectively zero error.
Table 8.1. Error rates for tgamma
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Sun compiler version 0x5150 
Microsoft Visual C++ version 14.1 


factorials 
Max = 0ε (Mean = 0ε) 
Max = 2.67ε (Mean = 0.617ε) 
Max = 172ε (Mean = 41ε) 
Max = 1.85ε (Mean = 0.566ε) 
near 0 
Max = 0ε (Mean = 0ε) 
Max = 2ε (Mean = 0.608ε) 
Max = 2ε (Mean = 0.647ε) 
Max = 1.5ε (Mean = 0.635ε) 
near 1 
Max = 0ε (Mean = 0ε) 
Max = 2.51ε (Mean = 1.02ε) 
Max = 3.01ε (Mean = 1.06ε) 
Max = 1.1ε (Mean = 0.59ε) 
near 2 
Max = 0ε (Mean = 0ε) 
Max = 4.1ε (Mean = 1.55ε) 
Max = 5.01ε (Mean = 1.89ε) 
Max = 2ε (Mean = 0.733ε) 
near 10 
Max = 0ε (Mean = 0ε) 
Max = 1.75ε (Mean = 0.895ε) 
Max = 1.75ε (Mean = 0.819ε) 
Max = 1.86ε (Mean = 0.881ε) 
near 55 
Max = 0ε (Mean = 0ε) 
Max = 2.69ε (Mean = 1.09ε) 
Max = 98.5ε (Mean = 53.4ε) 
Max = 2.7ε (Mean = 1.35ε) 
Table 8.2. Error rates for tgamma1pm1
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Sun compiler version 0x5150 
Microsoft Visual C++ version 14.1 


tgamma1pm1(dz) 
Max = 0ε (Mean = 0ε) 
Max = 1.12ε (Mean = 0.49ε) 
Max = 6.61ε (Mean = 0.84ε) 
Max = 3.31ε (Mean = 0.517ε) 
The following error plot are based on an exhaustive search of the functions
domain, MSVC15.5 at double
precision, and GCC7.1/Ubuntu for long
double
and __float128
.
The gamma is relatively easy to test: factorials and halfinteger 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.