boost_1_45_0/libs/math/doc/sf_and_dist/lgamma.qbk - nest-learning-thermostat/5.0/boost - Git at Google

 [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).
 ]
	[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).
	]