You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and dots ('.'), can be up to 35 characters long. Letters must be lowercase.
450 lines
14 KiB
450 lines
14 KiB
// (C) Copyright John Maddock 2006. |
|
// Use, modification and distribution are subject to 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) |
|
|
|
#ifndef BOOST_MATH_SF_DIGAMMA_HPP |
|
#define BOOST_MATH_SF_DIGAMMA_HPP |
|
|
|
#ifdef _MSC_VER |
|
#pragma once |
|
#endif |
|
|
|
#include <boost/math/tools/rational.hpp> |
|
#include <boost/math/tools/promotion.hpp> |
|
#include <boost/math/policies/error_handling.hpp> |
|
#include <boost/math/constants/constants.hpp> |
|
#include <boost/mpl/comparison.hpp> |
|
|
|
namespace boost{ |
|
namespace math{ |
|
namespace detail{ |
|
// |
|
// Begin by defining the smallest value for which it is safe to |
|
// use the asymptotic expansion for digamma: |
|
// |
|
inline unsigned digamma_large_lim(const mpl::int_<0>*) |
|
{ return 20; } |
|
|
|
inline unsigned digamma_large_lim(const void*) |
|
{ return 10; } |
|
// |
|
// Implementations of the asymptotic expansion come next, |
|
// the coefficients of the series have been evaluated |
|
// in advance at high precision, and the series truncated |
|
// at the first term that's too small to effect the result. |
|
// Note that the series becomes divergent after a while |
|
// so truncation is very important. |
|
// |
|
// This first one gives 34-digit precision for x >= 20: |
|
// |
|
template <class T> |
|
inline T digamma_imp_large(T x, const mpl::int_<0>*) |
|
{ |
|
BOOST_MATH_STD_USING // ADL of std functions. |
|
static const T P[] = { |
|
0.083333333333333333333333333333333333333333333333333L, |
|
-0.0083333333333333333333333333333333333333333333333333L, |
|
0.003968253968253968253968253968253968253968253968254L, |
|
-0.0041666666666666666666666666666666666666666666666667L, |
|
0.0075757575757575757575757575757575757575757575757576L, |
|
-0.021092796092796092796092796092796092796092796092796L, |
|
0.083333333333333333333333333333333333333333333333333L, |
|
-0.44325980392156862745098039215686274509803921568627L, |
|
3.0539543302701197438039543302701197438039543302701L, |
|
-26.456212121212121212121212121212121212121212121212L, |
|
281.4601449275362318840579710144927536231884057971L, |
|
-3607.510546398046398046398046398046398046398046398L, |
|
54827.583333333333333333333333333333333333333333333L, |
|
-974936.82385057471264367816091954022988505747126437L, |
|
20052695.796688078946143462272494530559046688078946L, |
|
-472384867.72162990196078431372549019607843137254902L, |
|
12635724795.916666666666666666666666666666666666667L |
|
}; |
|
x -= 1; |
|
T result = log(x); |
|
result += 1 / (2 * x); |
|
T z = 1 / (x*x); |
|
result -= z * tools::evaluate_polynomial(P, z); |
|
return result; |
|
} |
|
// |
|
// 19-digit precision for x >= 10: |
|
// |
|
template <class T> |
|
inline T digamma_imp_large(T x, const mpl::int_<64>*) |
|
{ |
|
BOOST_MATH_STD_USING // ADL of std functions. |
|
static const T P[] = { |
|
0.083333333333333333333333333333333333333333333333333L, |
|
-0.0083333333333333333333333333333333333333333333333333L, |
|
0.003968253968253968253968253968253968253968253968254L, |
|
-0.0041666666666666666666666666666666666666666666666667L, |
|
0.0075757575757575757575757575757575757575757575757576L, |
|
-0.021092796092796092796092796092796092796092796092796L, |
|
0.083333333333333333333333333333333333333333333333333L, |
|
-0.44325980392156862745098039215686274509803921568627L, |
|
3.0539543302701197438039543302701197438039543302701L, |
|
-26.456212121212121212121212121212121212121212121212L, |
|
281.4601449275362318840579710144927536231884057971L, |
|
}; |
|
x -= 1; |
|
T result = log(x); |
|
result += 1 / (2 * x); |
|
T z = 1 / (x*x); |
|
result -= z * tools::evaluate_polynomial(P, z); |
|
return result; |
|
} |
|
// |
|
// 17-digit precision for x >= 10: |
|
// |
|
template <class T> |
|
inline T digamma_imp_large(T x, const mpl::int_<53>*) |
|
{ |
|
BOOST_MATH_STD_USING // ADL of std functions. |
|
static const T P[] = { |
|
0.083333333333333333333333333333333333333333333333333L, |
|
-0.0083333333333333333333333333333333333333333333333333L, |
|
0.003968253968253968253968253968253968253968253968254L, |
|
-0.0041666666666666666666666666666666666666666666666667L, |
|
0.0075757575757575757575757575757575757575757575757576L, |
|
-0.021092796092796092796092796092796092796092796092796L, |
|
0.083333333333333333333333333333333333333333333333333L, |
|
-0.44325980392156862745098039215686274509803921568627L |
|
}; |
|
x -= 1; |
|
T result = log(x); |
|
result += 1 / (2 * x); |
|
T z = 1 / (x*x); |
|
result -= z * tools::evaluate_polynomial(P, z); |
|
return result; |
|
} |
|
// |
|
// 9-digit precision for x >= 10: |
|
// |
|
template <class T> |
|
inline T digamma_imp_large(T x, const mpl::int_<24>*) |
|
{ |
|
BOOST_MATH_STD_USING // ADL of std functions. |
|
static const T P[] = { |
|
0.083333333333333333333333333333333333333333333333333L, |
|
-0.0083333333333333333333333333333333333333333333333333L, |
|
0.003968253968253968253968253968253968253968253968254L |
|
}; |
|
x -= 1; |
|
T result = log(x); |
|
result += 1 / (2 * x); |
|
T z = 1 / (x*x); |
|
result -= z * tools::evaluate_polynomial(P, z); |
|
return result; |
|
} |
|
// |
|
// Now follow rational approximations over the range [1,2]. |
|
// |
|
// 35-digit precision: |
|
// |
|
template <class T> |
|
T digamma_imp_1_2(T x, const mpl::int_<0>*) |
|
{ |
|
// |
|
// Now the approximation, we use the form: |
|
// |
|
// digamma(x) = (x - root) * (Y + R(x-1)) |
|
// |
|
// Where root is the location of the positive root of digamma, |
|
// Y is a constant, and R is optimised for low absolute error |
|
// compared to Y. |
|
// |
|
// Max error found at 128-bit long double precision: 5.541e-35 |
|
// Maximum Deviation Found (approximation error): 1.965e-35 |
|
// |
|
static const float Y = 0.99558162689208984375F; |
|
|
|
static const T root1 = 1569415565.0 / 1073741824uL; |
|
static const T root2 = (381566830.0 / 1073741824uL) / 1073741824uL; |
|
static const T root3 = ((111616537.0 / 1073741824uL) / 1073741824uL) / 1073741824uL; |
|
static const T root4 = (((503992070.0 / 1073741824uL) / 1073741824uL) / 1073741824uL) / 1073741824uL; |
|
static const T root5 = 0.52112228569249997894452490385577338504019838794544e-36L; |
|
|
|
static const T P[] = { |
|
0.25479851061131551526977464225335883769L, |
|
-0.18684290534374944114622235683619897417L, |
|
-0.80360876047931768958995775910991929922L, |
|
-0.67227342794829064330498117008564270136L, |
|
-0.26569010991230617151285010695543858005L, |
|
-0.05775672694575986971640757748003553385L, |
|
-0.0071432147823164975485922555833274240665L, |
|
-0.00048740753910766168912364555706064993274L, |
|
-0.16454996865214115723416538844975174761e-4L, |
|
-0.20327832297631728077731148515093164955e-6L |
|
}; |
|
static const T Q[] = { |
|
1, |
|
2.6210924610812025425088411043163287646L, |
|
2.6850757078559596612621337395886392594L, |
|
1.4320913706209965531250495490639289418L, |
|
0.4410872083455009362557012239501953402L, |
|
0.081385727399251729505165509278152487225L, |
|
0.0089478633066857163432104815183858149496L, |
|
0.00055861622855066424871506755481997374154L, |
|
0.1760168552357342401304462967950178554e-4L, |
|
0.20585454493572473724556649516040874384e-6L, |
|
-0.90745971844439990284514121823069162795e-11L, |
|
0.48857673606545846774761343500033283272e-13L, |
|
}; |
|
T g = x - root1; |
|
g -= root2; |
|
g -= root3; |
|
g -= root4; |
|
g -= root5; |
|
T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1)); |
|
T result = g * Y + g * r; |
|
|
|
return result; |
|
} |
|
// |
|
// 19-digit precision: |
|
// |
|
template <class T> |
|
T digamma_imp_1_2(T x, const mpl::int_<64>*) |
|
{ |
|
// |
|
// Now the approximation, we use the form: |
|
// |
|
// digamma(x) = (x - root) * (Y + R(x-1)) |
|
// |
|
// Where root is the location of the positive root of digamma, |
|
// Y is a constant, and R is optimised for low absolute error |
|
// compared to Y. |
|
// |
|
// Max error found at 80-bit long double precision: 5.016e-20 |
|
// Maximum Deviation Found (approximation error): 3.575e-20 |
|
// |
|
static const float Y = 0.99558162689208984375F; |
|
|
|
static const T root1 = 1569415565.0 / 1073741824uL; |
|
static const T root2 = (381566830.0 / 1073741824uL) / 1073741824uL; |
|
static const T root3 = 0.9016312093258695918615325266959189453125e-19L; |
|
|
|
static const T P[] = { |
|
0.254798510611315515235L, |
|
-0.314628554532916496608L, |
|
-0.665836341559876230295L, |
|
-0.314767657147375752913L, |
|
-0.0541156266153505273939L, |
|
-0.00289268368333918761452L |
|
}; |
|
static const T Q[] = { |
|
1, |
|
2.1195759927055347547L, |
|
1.54350554664961128724L, |
|
0.486986018231042975162L, |
|
0.0660481487173569812846L, |
|
0.00298999662592323990972L, |
|
-0.165079794012604905639e-5L, |
|
0.317940243105952177571e-7L |
|
}; |
|
T g = x - root1; |
|
g -= root2; |
|
g -= root3; |
|
T r = tools::evaluate_polynomial(P, x-1) / tools::evaluate_polynomial(Q, x-1); |
|
T result = g * Y + g * r; |
|
|
|
return result; |
|
} |
|
// |
|
// 18-digit precision: |
|
// |
|
template <class T> |
|
T digamma_imp_1_2(T x, const mpl::int_<53>*) |
|
{ |
|
// |
|
// Now the approximation, we use the form: |
|
// |
|
// digamma(x) = (x - root) * (Y + R(x-1)) |
|
// |
|
// Where root is the location of the positive root of digamma, |
|
// Y is a constant, and R is optimised for low absolute error |
|
// compared to Y. |
|
// |
|
// Maximum Deviation Found: 1.466e-18 |
|
// At double precision, max error found: 2.452e-17 |
|
// |
|
static const float Y = 0.99558162689208984F; |
|
|
|
static const T root1 = 1569415565.0 / 1073741824uL; |
|
static const T root2 = (381566830.0 / 1073741824uL) / 1073741824uL; |
|
static const T root3 = 0.9016312093258695918615325266959189453125e-19L; |
|
|
|
static const T P[] = { |
|
0.25479851061131551L, |
|
-0.32555031186804491L, |
|
-0.65031853770896507L, |
|
-0.28919126444774784L, |
|
-0.045251321448739056L, |
|
-0.0020713321167745952L |
|
}; |
|
static const T Q[] = { |
|
1L, |
|
2.0767117023730469L, |
|
1.4606242909763515L, |
|
0.43593529692665969L, |
|
0.054151797245674225L, |
|
0.0021284987017821144L, |
|
-0.55789841321675513e-6L |
|
}; |
|
T g = x - root1; |
|
g -= root2; |
|
g -= root3; |
|
T r = tools::evaluate_polynomial(P, x-1) / tools::evaluate_polynomial(Q, x-1); |
|
T result = g * Y + g * r; |
|
|
|
return result; |
|
} |
|
// |
|
// 9-digit precision: |
|
// |
|
template <class T> |
|
inline T digamma_imp_1_2(T x, const mpl::int_<24>*) |
|
{ |
|
// |
|
// Now the approximation, we use the form: |
|
// |
|
// digamma(x) = (x - root) * (Y + R(x-1)) |
|
// |
|
// Where root is the location of the positive root of digamma, |
|
// Y is a constant, and R is optimised for low absolute error |
|
// compared to Y. |
|
// |
|
// Maximum Deviation Found: 3.388e-010 |
|
// At float precision, max error found: 2.008725e-008 |
|
// |
|
static const float Y = 0.99558162689208984f; |
|
static const T root = 1532632.0f / 1048576; |
|
static const T root_minor = static_cast<T>(0.3700660185912626595423257213284682051735604e-6L); |
|
static const T P[] = { |
|
0.25479851023250261e0, |
|
-0.44981331915268368e0, |
|
-0.43916936919946835e0, |
|
-0.61041765350579073e-1 |
|
}; |
|
static const T Q[] = { |
|
0.1e1, |
|
0.15890202430554952e1, |
|
0.65341249856146947e0, |
|
0.63851690523355715e-1 |
|
}; |
|
T g = x - root; |
|
g -= root_minor; |
|
T r = tools::evaluate_polynomial(P, x-1) / tools::evaluate_polynomial(Q, x-1); |
|
T result = g * Y + g * r; |
|
|
|
return result; |
|
} |
|
|
|
template <class T, class Tag, class Policy> |
|
T digamma_imp(T x, const Tag* t, const Policy& pol) |
|
{ |
|
// |
|
// This handles reflection of negative arguments, and all our |
|
// error handling, then forwards to the T-specific approximation. |
|
// |
|
BOOST_MATH_STD_USING // ADL of std functions. |
|
|
|
T result = 0; |
|
// |
|
// Check for negative arguments and use reflection: |
|
// |
|
if(x < 0) |
|
{ |
|
// Reflect: |
|
x = 1 - x; |
|
// Argument reduction for tan: |
|
T remainder = x - floor(x); |
|
// Shift to negative if > 0.5: |
|
if(remainder > 0.5) |
|
{ |
|
remainder -= 1; |
|
} |
|
// |
|
// check for evaluation at a negative pole: |
|
// |
|
if(remainder == 0) |
|
{ |
|
return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, (1-x), pol); |
|
} |
|
result = constants::pi<T>() / tan(constants::pi<T>() * remainder); |
|
} |
|
// |
|
// If we're above the lower-limit for the |
|
// asymptotic expansion then use it: |
|
// |
|
if(x >= digamma_large_lim(t)) |
|
{ |
|
result += digamma_imp_large(x, t); |
|
} |
|
else |
|
{ |
|
// |
|
// If x > 2 reduce to the interval [1,2]: |
|
// |
|
while(x > 2) |
|
{ |
|
x -= 1; |
|
result += 1/x; |
|
} |
|
// |
|
// If x < 1 use recurrance to shift to > 1: |
|
// |
|
if(x < 1) |
|
{ |
|
result = -1/x; |
|
x += 1; |
|
} |
|
result += digamma_imp_1_2(x, t); |
|
} |
|
return result; |
|
} |
|
|
|
} // namespace detail |
|
|
|
template <class T, class Policy> |
|
inline typename tools::promote_args<T>::type |
|
digamma(T x, const Policy& pol) |
|
{ |
|
typedef typename tools::promote_args<T>::type result_type; |
|
typedef typename policies::evaluation<result_type, Policy>::type value_type; |
|
typedef typename policies::precision<T, Policy>::type precision_type; |
|
typedef typename mpl::if_< |
|
mpl::or_< |
|
mpl::less_equal<precision_type, mpl::int_<0> >, |
|
mpl::greater<precision_type, mpl::int_<64> > |
|
>, |
|
mpl::int_<0>, |
|
typename mpl::if_< |
|
mpl::less<precision_type, mpl::int_<25> >, |
|
mpl::int_<24>, |
|
typename mpl::if_< |
|
mpl::less<precision_type, mpl::int_<54> >, |
|
mpl::int_<53>, |
|
mpl::int_<64> |
|
>::type |
|
>::type |
|
>::type tag_type; |
|
|
|
return policies::checked_narrowing_cast<result_type, Policy>(detail::digamma_imp( |
|
static_cast<value_type>(x), |
|
static_cast<const tag_type*>(0), pol), "boost::math::digamma<%1%>(%1%)"); |
|
} |
|
|
|
template <class T> |
|
inline typename tools::promote_args<T>::type |
|
digamma(T x) |
|
{ |
|
return digamma(x, policies::policy<>()); |
|
} |
|
|
|
} // namespace math |
|
} // namespace boost |
|
#endif |
|
|
|
|