| [section:lgamma Log Gamma] |
| |
| [h4 Synopsis] |
| |
| `` |
| #include <boost/math/special_functions/gamma.hpp> |
| `` |
| |
| namespace boost{ namespace math{ |
| |
| template <class T> |
| ``__sf_result`` lgamma(T z); |
| |
| template <class T, class ``__Policy``> |
| ``__sf_result`` lgamma(T z, const ``__Policy``&); |
| |
| template <class T> |
| ``__sf_result`` lgamma(T z, int* sign); |
| |
| template <class T, class ``__Policy``> |
| ``__sf_result`` lgamma(T z, int* sign, const ``__Policy``&); |
| |
| }} // namespaces |
| |
| [h4 Description] |
| |
| The [@http://en.wikipedia.org/wiki/Gamma_function lgamma function] is defined by: |
| |
| [equation lgamm1] |
| |
| The second form of the function takes a pointer to an integer, |
| which if non-null is set on output to the sign of tgamma(z). |
| |
| [optional_policy] |
| |
| [graph lgamma] |
| |
| There are effectively two versions of this 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. In practice, any built-in |
| floating-point type you will encounter has an appropriate __lanczos |
| defined for it. It is also possible, given enough machine time, to generate |
| further __lanczos's using the program libs/math/tools/lanczos_generator.cpp. |
| |
| The return type of these functions is computed using the __arg_pomotion_rules: |
| the result is of type `double` if T is an integer type, or type T otherwise. |
| |
| [h4 Accuracy] |
| |
| 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, __glibc, __hpc and |
| __cephes libraries. Unless otherwise specified any |
| floating point type that is narrower than the one shown will have |
| __zero_error. |
| |
| Note that while the relative errors near the positive roots of lgamma |
| are very low, the lgamma function has an infinite number of irrational |
| roots for negative arguments: very close to these negative roots only |
| a low absolute error can be guaranteed. |
| |
| [table |
| [[Significand Size] [Platform and Compiler] [Factorials and Half factorials] [Values Near Zero] [Values Near 1 or 2] [Values Near a Negative Pole]] |
| [[53] [Win32 Visual C++ 8] |
| [Peak=0.88 Mean=0.14 |
| |
| (GSL=33) (__cephes=1.5)] |
| [Peak=0.96 Mean=0.46 |
| |
| (GSL=5.2) (__cephes=1.1)] |
| [Peak=0.86 Mean=0.46 |
| |
| (GSL=1168) (__cephes~500000)] |
| [Peak=4.2 Mean=1.3 |
| |
| (GSL=25) (__cephes=1.6)] ] |
| [[64] [Linux IA32 / GCC] |
| [Peak=1.9 Mean=0.43 |
| |
| (__glibc Peak=1.7 Mean=0.49)] |
| [Peak=1.4 Mean=0.57 |
| |
| (__glibc Peak= 0.96 Mean=0.54)] |
| [Peak=0.86 Mean=0.35 |
| |
| (__glibc Peak=0.74 Mean=0.26)] |
| |
| [Peak=6.0 Mean=1.8 |
| |
| (__glibc Peak=3.0 Mean=0.86)] ] |
| [[64] [Linux IA64 / GCC] |
| [Peak=0.99 Mean=0.12 |
| |
| (__glibc Peak 0)] |
| |
| [Pek=1.2 Mean=0.6 |
| |
| (__glibc Peak 0)] |
| [Peak=0.86 Mean=0.16 |
| |
| (__glibc Peak 0)] |
| [Peak=2.3 Mean=0.69 |
| |
| (__glibc Peak 0)] ] |
| [[113] [HPUX IA64, aCC A.06.06] |
| [Peak=0.96 Mean=0.13 |
| |
| (__hpc Peak 0)] |
| [Peak=0.99 Mean=0.53 |
| |
| (__hpc Peak 0)] |
| [Peak=0.9 Mean=0.4 |
| |
| (__hpc Peak 0)] |
| [Peak=3.0 Mean=0.9 |
| |
| (__hpc Peak 0)] ] |
| ] |
| |
| [h4 Testing] |
| |
| The main tests for this function involve comparisons against the logs of |
| the factorials which can be independently calculated to very high accuracy. |
| |
| Random tests in key problem areas are also used. |
| |
| [h4 Implementation] |
| |
| The generic version of this function is implemented by combining the series and |
| continued fraction representations for the incomplete gamma function: |
| |
| [equation lgamm2] |
| |
| where /l/ is an arbitrary integration limit: choosing [^l = max(10, a)] |
| seems to work fairly well. For negative /z/ the logarithm version of the |
| reflection formula is used: |
| |
| [equation lgamm3] |
| |
| For types of known precision, the __lanczos is used, a traits class |
| `boost::math::lanczos::lanczos_traits` maps type T to an appropriate |
| approximation. The logarithmic version of the __lanczos is: |
| |
| [equation lgamm4] |
| |
| Where L[sub e,g][space] is the Lanczos sum, scaled by e[super g]. |
| |
| As before the reflection formula is used for /z < 0/. |
| |
| When z is very near 1 or 2, then the logarithmic version of the __lanczos |
| suffers very badly from cancellation error: indeed for values sufficiently |
| close to 1 or 2, arbitrarily large relative errors can be obtained (even though |
| the absolute error is tiny). |
| |
| For types with up to 113 bits of precision |
| (up to and including 128-bit long doubles), root-preserving |
| rational approximations [jm_rationals] are used |
| over the intervals [1,2] and [2,3]. Over the interval [2,3] the approximation |
| form used is: |
| |
| lgamma(z) = (z-2)(z+1)(Y + R(z-2)); |
| |
| Where Y is a constant, and R(z-2) is the rational approximation: optimised |
| so that it's absolute error is tiny compared to Y. In addition |
| small values of z greater |
| than 3 can handled by argument reduction using the recurrence relation: |
| |
| lgamma(z+1) = log(z) + lgamma(z); |
| |
| Over the interval [1,2] two approximations have to be used, one for small z uses: |
| |
| lgamma(z) = (z-1)(z-2)(Y + R(z-1)); |
| |
| Once again Y is a constant, and R(z-1) is optimised for low absolute error |
| compared to Y. For z > 1.5 the above form wouldn't converge to a |
| minimax solution but this similar form does: |
| |
| lgamma(z) = (2-z)(1-z)(Y + R(2-z)); |
| |
| Finally for z < 1 the recurrence relation can be used to move to z > 1: |
| |
| lgamma(z) = lgamma(z+1) - log(z); |
| |
| Note that while this involves a subtraction, it appears not |
| to suffer from cancellation error: as z decreases from 1 |
| the `-log(z)` term grows positive much more |
| rapidly than the `lgamma(z+1)` term becomes negative. So in this |
| specific case, significant digits are preserved, rather than cancelled. |
| |
| For other types which do have a __lanczos defined for them |
| the current solution is as follows: imagine we |
| balance the two terms in the __lanczos by dividing the power term by its value |
| at /z = 1/, and then multiplying the Lanczos coefficients by the same value. |
| Now each term will take the value 1 at /z = 1/ and we can rearrange the power terms |
| in terms of log1p. Likewise if we subtract 1 from the Lanczos sum part |
| (algebraically, by subtracting the value of each term at /z = 1/), we obtain |
| a new summation that can be also be fed into log1p. Crucially, all of the |
| terms tend to zero, as /z -> 1/: |
| |
| [equation lgamm5] |
| |
| The C[sub k][space] terms in the above are the same as in the __lanczos. |
| |
| A similar rearrangement can be performed at /z = 2/: |
| |
| [equation lgamm6] |
| |
| [endsect][/section:lgamma The Log Gamma Function] |
| |
| [/ |
| Copyright 2006 John Maddock and Paul A. Bristow. |
| Distributed under the Boost Software License, Version 1.0. |
| (See accompanying file LICENSE_1_0.txt or copy at |
| http://www.boost.org/LICENSE_1_0.txt). |
| ] |