mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 04:22:09 +00:00
Rewritten tgamma1pm1 in terms of lgamma rational approximations: now more efficient and more accurate!
Minor patch to exponential distribution tests. Removed dead code from Remez app. [SVN r3331]
This commit is contained in:
John Maddock
@@ -175,7 +175,7 @@ T gamma_imp(T z, const L& l)
|
||||
// lgamma for small arguments:
|
||||
//
|
||||
template <class T, class L>
|
||||
T lgamma_small_imp(T z, const mpl::int_<64>&, const L& /* l */)
|
||||
T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<64>&, const L& /* l */)
|
||||
{
|
||||
// This version uses rational approximations for small
|
||||
// values of z accurate enough for 64-bit mantissas
|
||||
@@ -188,7 +188,7 @@ T lgamma_small_imp(T z, const mpl::int_<64>&, const L& /* l */)
|
||||
{
|
||||
result = -log(z);
|
||||
}
|
||||
else if((z == 1) || (z == 2))
|
||||
else if((zm1 == 0) || (zm2 == 0))
|
||||
{
|
||||
// nothing to do, result is zero....
|
||||
}
|
||||
@@ -198,10 +198,16 @@ T lgamma_small_imp(T z, const mpl::int_<64>&, const L& /* l */)
|
||||
// Begin by performing argument reduction until
|
||||
// z is in [2,3):
|
||||
//
|
||||
while(z >= 3)
|
||||
if(z >= 3)
|
||||
{
|
||||
z -= 1;
|
||||
result += log(z);
|
||||
do
|
||||
{
|
||||
z -= 1;
|
||||
zm2 -= 1;
|
||||
result += log(z);
|
||||
}while(z >= 3);
|
||||
// Update zm2, we need it below:
|
||||
zm2 = z - 2;
|
||||
}
|
||||
|
||||
//
|
||||
@@ -244,7 +250,6 @@ T lgamma_small_imp(T z, const mpl::int_<64>&, const L& /* l */)
|
||||
|
||||
static const float Y = 0.158963680267333984375f;
|
||||
|
||||
T zm2 = z - 2;
|
||||
T r = zm2 * (z + 1);
|
||||
T R = tools::evaluate_polynomial(P, zm2);
|
||||
R /= tools::evaluate_polynomial(Q, zm2);
|
||||
@@ -257,15 +262,13 @@ T lgamma_small_imp(T z, const mpl::int_<64>&, const L& /* l */)
|
||||
// If z is less than 1 use recurrance to shift to
|
||||
// z in the interval [1,2]:
|
||||
//
|
||||
T zm1;
|
||||
if(z < 1)
|
||||
{
|
||||
result += -log(z);
|
||||
zm2 = zm1;
|
||||
zm1 = z;
|
||||
z += 1;
|
||||
}
|
||||
else
|
||||
zm1 = z - 1;
|
||||
//
|
||||
// Two approximations, on for z in [1,1.5] and
|
||||
// one for z in [1.5,2]:
|
||||
@@ -309,7 +312,7 @@ T lgamma_small_imp(T z, const mpl::int_<64>&, const L& /* l */)
|
||||
};
|
||||
|
||||
T r = tools::evaluate_polynomial(P, zm1) / tools::evaluate_polynomial(Q, zm1);
|
||||
T prefix = zm1 * (z - 2);
|
||||
T prefix = zm1 * zm2;
|
||||
|
||||
result += prefix * Y + prefix * r;
|
||||
}
|
||||
@@ -349,10 +352,8 @@ T lgamma_small_imp(T z, const mpl::int_<64>&, const L& /* l */)
|
||||
-0.00100122979330632606976L,
|
||||
-0.818096163862371810295e-6L
|
||||
};
|
||||
// (2 - x) * (1 - x) * (c + R(2 - x))
|
||||
T zm2 = 2 - z;
|
||||
T r = -zm2 * zm1;
|
||||
T R = tools::evaluate_polynomial(P, zm2) / tools::evaluate_polynomial(Q, zm2);
|
||||
T r = zm2 * zm1;
|
||||
T R = tools::evaluate_polynomial(P, -zm2) / tools::evaluate_polynomial(Q, -zm2);
|
||||
|
||||
result += r * Y + r * R;
|
||||
}
|
||||
@@ -360,7 +361,7 @@ T lgamma_small_imp(T z, const mpl::int_<64>&, const L& /* l */)
|
||||
return result;
|
||||
}
|
||||
template <class T, class L>
|
||||
T lgamma_small_imp(T z, const mpl::int_<113>&, const L& /* l */)
|
||||
T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<113>&, const L& /* l */)
|
||||
{
|
||||
//
|
||||
// This version uses rational approximations for small
|
||||
@@ -373,7 +374,7 @@ T lgamma_small_imp(T z, const mpl::int_<113>&, const L& /* l */)
|
||||
{
|
||||
result = -log(z);
|
||||
}
|
||||
else if((z == 1) || (z == 2))
|
||||
else if((zm1 == 0) || (zm2 == 0))
|
||||
{
|
||||
// nothing to do, result is zero....
|
||||
}
|
||||
@@ -383,10 +384,14 @@ T lgamma_small_imp(T z, const mpl::int_<113>&, const L& /* l */)
|
||||
// Begin by performing argument reduction until
|
||||
// z is in [2,3):
|
||||
//
|
||||
while(z >= 3)
|
||||
if(z >= 3)
|
||||
{
|
||||
z -= 1;
|
||||
result += log(z);
|
||||
do
|
||||
{
|
||||
z -= 1;
|
||||
result += log(z);
|
||||
}while(z >= 3);
|
||||
zm2 = z - 2;
|
||||
}
|
||||
|
||||
//
|
||||
@@ -399,7 +404,6 @@ T lgamma_small_imp(T z, const mpl::int_<113>&, const L& /* l */)
|
||||
// is small compared to the constant Y - then any rounding
|
||||
// error in it's computation will get wiped out.
|
||||
//
|
||||
T zm2 = z - 2;
|
||||
if(zm2 <= 0.5)
|
||||
{
|
||||
// Max error found at 128-bit long double precision 2.079e-36
|
||||
@@ -472,7 +476,7 @@ T lgamma_small_imp(T z, const mpl::int_<113>&, const L& /* l */)
|
||||
-0.21417746406796251935895373591106602902e-7L
|
||||
};
|
||||
|
||||
T R = tools::evaluate_polynomial(P, z - 2.5) / tools::evaluate_polynomial(Q, z - 2.5);
|
||||
T R = tools::evaluate_polynomial(P, zm2 - 0.5) / tools::evaluate_polynomial(Q, zm2 - 0.5);
|
||||
|
||||
static const float Y = 0.16830921173095703125f;
|
||||
|
||||
@@ -487,15 +491,13 @@ T lgamma_small_imp(T z, const mpl::int_<113>&, const L& /* l */)
|
||||
// If z is less than 1 use recurrance to shift to
|
||||
// z in the interval [1,2]:
|
||||
//
|
||||
T zm1;
|
||||
if(z < 1)
|
||||
{
|
||||
result += -log(z);
|
||||
zm2 = zm1;
|
||||
zm1 = z;
|
||||
z += 1;
|
||||
}
|
||||
else
|
||||
zm1 = z - 1;
|
||||
//
|
||||
// Three approximations, on for z in [1,1.35], [1.35,1.625] and [1.625,1]
|
||||
//
|
||||
@@ -547,7 +549,7 @@ T lgamma_small_imp(T z, const mpl::int_<113>&, const L& /* l */)
|
||||
};
|
||||
|
||||
T r = tools::evaluate_polynomial(P, zm1) / tools::evaluate_polynomial(Q, zm1);
|
||||
T prefix = zm1 * (z - 2);
|
||||
T prefix = zm1 * zm2;
|
||||
|
||||
result += prefix * Y + prefix * r;
|
||||
}
|
||||
@@ -595,10 +597,8 @@ T lgamma_small_imp(T z, const mpl::int_<113>&, const L& /* l */)
|
||||
0.00013288854760548251757651556792598235735L,
|
||||
-0.17194794958274081373243161848194745111e-5L
|
||||
};
|
||||
// (2 - x) * (1 - x) * (c + R(2 - x))
|
||||
T zm2 = 2 - z;
|
||||
T r = -zm2 * zm1;
|
||||
T R = tools::evaluate_polynomial(P, 1.625 - z) / tools::evaluate_polynomial(Q, 1.625 - z);
|
||||
T r = zm2 * zm1;
|
||||
T R = tools::evaluate_polynomial(P, 0.625 - zm1) / tools::evaluate_polynomial(Q, 0.625 - zm1);
|
||||
|
||||
result += r * Y + r * R;
|
||||
}
|
||||
@@ -638,9 +638,8 @@ T lgamma_small_imp(T z, const mpl::int_<113>&, const L& /* l */)
|
||||
-0.30202973883316730694433702165188835331e-6L
|
||||
};
|
||||
// (2 - x) * (1 - x) * (c + R(2 - x))
|
||||
T zm2 = 2 - z;
|
||||
T r = -zm2 * zm1;
|
||||
T R = tools::evaluate_polynomial(P, zm2) / tools::evaluate_polynomial(Q, zm2);
|
||||
T r = zm2 * zm1;
|
||||
T R = tools::evaluate_polynomial(P, -zm2) / tools::evaluate_polynomial(Q, -zm2);
|
||||
|
||||
result += r * Y + r * R;
|
||||
}
|
||||
@@ -648,7 +647,7 @@ T lgamma_small_imp(T z, const mpl::int_<113>&, const L& /* l */)
|
||||
return result;
|
||||
}
|
||||
template <class T, class L>
|
||||
T lgamma_small_imp(T z, const mpl::int_<0>&, const L& l)
|
||||
T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<0>&, const L& l)
|
||||
{
|
||||
//
|
||||
// No rational approximations are available because either
|
||||
@@ -676,15 +675,15 @@ T lgamma_small_imp(T z, const mpl::int_<0>&, const L& l)
|
||||
else if(z >= 1.5)
|
||||
{
|
||||
// special case near 2:
|
||||
T dz = z - 2;
|
||||
T dz = zm2;
|
||||
result = dz * log((z + L::g() - T(0.5)) / boost::math::constants::e<T>());
|
||||
result += boost::math::log1p(dz / (L::g() + T(1.5))) * T(1.5);
|
||||
result += boost::math::log1p(L::lanczos_sum_near_2(z));
|
||||
result += boost::math::log1p(L::lanczos_sum_near_2(dz));
|
||||
}
|
||||
else
|
||||
{
|
||||
// special case near 1:
|
||||
T dz = z - 1;
|
||||
T dz = zm1;
|
||||
result = dz * log((z + L::g() - T(0.5)) / boost::math::constants::e<T>());
|
||||
result += boost::math::log1p(dz / (L::g() + T(0.5))) / 2;
|
||||
result += boost::math::log1p(L::lanczos_sum_near_1(dz));
|
||||
@@ -737,7 +736,7 @@ T lgamma_imp(T z, const L& l, int* sign = 0)
|
||||
(std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::digits <= 113)),
|
||||
mpl::int_<113>, mpl::int_<0> >::type
|
||||
>::type tag_type;
|
||||
result = lgamma_small_imp(z, tag_type(), l);
|
||||
result = lgamma_small_imp(z, z - 1, z - 2, tag_type(), l);
|
||||
}
|
||||
else if((z >= 3) && (z < 100))
|
||||
{
|
||||
@@ -918,45 +917,80 @@ T lgamma_imp(T z, const lanczos::undefined_lanczos&, int*sign)
|
||||
// This helper calculates tgamma(dz+1)-1 without cancellation errors,
|
||||
// used by the upper incomplete gamma with z < 1:
|
||||
//
|
||||
template <class T, class L>
|
||||
T tgammap1m1_imp(T dz, const L&)
|
||||
template <class T, class Tag, class L>
|
||||
T tgammap1m1_imp(T dz, Tag const& tag, const L& l)
|
||||
{
|
||||
//
|
||||
using namespace std;
|
||||
|
||||
T zgh = (L::g() + T(0.5) + dz) / boost::math::constants::e<T>();
|
||||
T A = boost::math::powm1(zgh, dz);
|
||||
T B = boost::math::sqrt1pm1(dz / (L::g() + T(0.5)));
|
||||
T C = L::lanczos_sum_near_1(dz);
|
||||
T Ap1 = pow(zgh, dz);
|
||||
T Bp1 = sqrt(1 + (dz / (L::g() + T(0.5))) );
|
||||
T result;
|
||||
if(dz < 0)
|
||||
{
|
||||
if(dz < -0.5)
|
||||
{
|
||||
// Best method is simply to subtract 1 from tgamma:
|
||||
result = boost::math::tgamma(1+dz) - 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Use expm1 on lgamma:
|
||||
result = boost::math::expm1(-boost::math::log1p(dz)
|
||||
+ lgamma_small_imp(dz+2, dz + 1, dz, tag, l));
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if(dz < 2)
|
||||
{
|
||||
// Use expm1 on lgamma:
|
||||
result = boost::math::expm1(lgamma_small_imp(dz+1, dz, dz-1, tag, l));
|
||||
}
|
||||
else
|
||||
{
|
||||
// Best method is simply to subtract 1 from tgamma:
|
||||
result = boost::math::tgamma(1+dz) - 1;
|
||||
}
|
||||
}
|
||||
|
||||
return Bp1 * (A + C * Ap1) + B;
|
||||
return result;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
T tgammap1m1_imp(T dz,
|
||||
template <class T, class Tag>
|
||||
T tgammap1m1_imp(T dz, Tag const& tag,
|
||||
const ::boost::math::lanczos::undefined_lanczos& l)
|
||||
{
|
||||
//
|
||||
// TODO FIXME!!!!
|
||||
// There should be a better solution than this, but the
|
||||
// algebra isn't easy for the general case....
|
||||
// Start by subracting 1 from tgamma:
|
||||
//
|
||||
T result = gamma_imp(1 + dz, l) - 1;
|
||||
if(std::pow(2.0, boost::math::tools::digits<T>()) * result
|
||||
< std::pow(2.0, std::numeric_limits<long double>::digits))
|
||||
//
|
||||
// Test the level of cancellation error observed: we loose one bit
|
||||
// for each power of 2 the result is less than 1. If we would get
|
||||
// more bits from our most precise lgamma rational approximation,
|
||||
// then use that instead:
|
||||
//
|
||||
if((dz > -0.5) && (dz < 2) && (ldexp(1.0, boost::math::tools::digits<T>()) * fabs(result) < 1e34))
|
||||
{
|
||||
// Cancellation errors mean that we have
|
||||
// fewer good digits left in the result than
|
||||
// there are digits in a long double. To limit
|
||||
// the damage, call the long double version:
|
||||
typedef boost::math::lanczos::lanczos_traits<long double>::evaluation_type eval_type;
|
||||
result = tgammap1m1_imp(boost::math::tools::real_cast<long double>(dz), eval_type());
|
||||
result = tgammap1m1_imp(dz, mpl::int_<113>(), boost::math::lanczos::lanczos24m113());
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
template <class T, class L>
|
||||
inline T tgammap1m1_imp(T dz, const L& l)
|
||||
{
|
||||
// Simple forwarding function.
|
||||
typedef typename mpl::if_c<
|
||||
(std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::digits <= 64)),
|
||||
mpl::int_<64>,
|
||||
typename mpl::if_c<
|
||||
(std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::digits <= 113)),
|
||||
mpl::int_<113>, mpl::int_<0> >::type
|
||||
>::type tag_type;
|
||||
return tgammap1m1_imp(dz, tag_type(), l);
|
||||
}
|
||||
|
||||
//
|
||||
// Series representation for upper fraction when z is small:
|
||||
//
|
||||
@@ -1609,7 +1643,15 @@ inline T tgamma1pm1(T z)
|
||||
BOOST_FPU_EXCEPTION_GUARD
|
||||
typedef typename lanczos::lanczos_traits<typename remove_cv<T>::type>::value_type value_type;
|
||||
typedef typename lanczos::lanczos_traits<typename remove_cv<T>::type>::evaluation_type evaluation_type;
|
||||
return tools::checked_narrowing_cast<typename remove_cv<T>::type>(detail::tgammap1m1_imp(static_cast<value_type>(z), evaluation_type()), BOOST_CURRENT_FUNCTION);
|
||||
typedef typename mpl::if_c<
|
||||
(std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::digits <= 64)),
|
||||
mpl::int_<64>,
|
||||
typename mpl::if_c<
|
||||
(std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::digits <= 113)),
|
||||
mpl::int_<113>, mpl::int_<0> >::type
|
||||
>::type tag_type;
|
||||
|
||||
return tools::checked_narrowing_cast<typename remove_cv<T>::type>(detail::tgammap1m1_imp(static_cast<value_type>(z), tag_type(), evaluation_type()), BOOST_CURRENT_FUNCTION);
|
||||
}
|
||||
|
||||
//
|
||||
|
||||
@@ -98,7 +98,7 @@ struct lanczos6
|
||||
}
|
||||
|
||||
template<class T>
|
||||
static T lanczos_sum_near_2(const T& z)
|
||||
static T lanczos_sum_near_2(const T& dz)
|
||||
{
|
||||
static const T d[5] = {
|
||||
static_cast<T>(5.748142489536043490764289256167080091892L),
|
||||
@@ -108,7 +108,7 @@ struct lanczos6
|
||||
static_cast<T>(0.002763134585812698552178368447708846850353L),
|
||||
};
|
||||
T result = 0;
|
||||
T dz = z - 2;
|
||||
T z = dz + 2;
|
||||
for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
|
||||
{
|
||||
result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
|
||||
@@ -215,7 +215,7 @@ struct lanczos11
|
||||
}
|
||||
|
||||
template<class T>
|
||||
static T lanczos_sum_near_2(const T& z)
|
||||
static T lanczos_sum_near_2(const T& dz)
|
||||
{
|
||||
static const T d[10] = {
|
||||
static_cast<T>(19.05889633808148715159575716844556056056L),
|
||||
@@ -230,7 +230,7 @@ struct lanczos11
|
||||
static_cast<T>(-0.493038376656195010308610694048822561263e-7L),
|
||||
};
|
||||
T result = 0;
|
||||
T dz = z - 2;
|
||||
T z = dz + 2;
|
||||
for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
|
||||
{
|
||||
result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
|
||||
@@ -347,7 +347,7 @@ struct lanczos13
|
||||
}
|
||||
|
||||
template<class T>
|
||||
static T lanczos_sum_near_2(const T& z)
|
||||
static T lanczos_sum_near_2(const T& dz)
|
||||
{
|
||||
static const T d[12] = {
|
||||
static_cast<T>(26.96979819614830698367887026728396466395L),
|
||||
@@ -364,7 +364,7 @@ struct lanczos13
|
||||
static_cast<T>(-0.9685385411006641478305219367315965391289e-9L),
|
||||
};
|
||||
T result = 0;
|
||||
T dz = z - 2;
|
||||
T z = z = 2;
|
||||
for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
|
||||
{
|
||||
result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
|
||||
@@ -526,7 +526,7 @@ struct lanczos22
|
||||
}
|
||||
|
||||
template<class T>
|
||||
static T lanczos_sum_near_2(const T& z)
|
||||
static T lanczos_sum_near_2(const T& dz)
|
||||
{
|
||||
static const T d[21] = {
|
||||
static_cast<T>(75.39272007105208086018421070699575462226L),
|
||||
@@ -552,7 +552,7 @@ struct lanczos22
|
||||
static_cast<T>(0.6580207998808093935798753964580596673177e-19L),
|
||||
};
|
||||
T result = 0;
|
||||
T dz = z - 2;
|
||||
T z = dz + 2;
|
||||
for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
|
||||
{
|
||||
result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
|
||||
@@ -634,7 +634,7 @@ struct lanczos6m24
|
||||
}
|
||||
|
||||
template<class T>
|
||||
static T lanczos_sum_near_2(const T& z)
|
||||
static T lanczos_sum_near_2(const T& dz)
|
||||
{
|
||||
static const T d[5] = {
|
||||
static_cast<T>(0.6534966888520080645505805298901130485464L),
|
||||
@@ -644,7 +644,7 @@ struct lanczos6m24
|
||||
static_cast<T>(-0.000750517993690428370380996157470900204524L),
|
||||
};
|
||||
T result = 0;
|
||||
T dz = z - 2;
|
||||
T z = dz + 2;
|
||||
for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
|
||||
{
|
||||
result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
|
||||
@@ -761,7 +761,7 @@ struct lanczos13m53
|
||||
}
|
||||
|
||||
template<class T>
|
||||
static T lanczos_sum_near_2(const T& z)
|
||||
static T lanczos_sum_near_2(const T& dz)
|
||||
{
|
||||
static const T d[12] = {
|
||||
static_cast<T>(6.565936202082889535528455955485877361223L),
|
||||
@@ -778,7 +778,7 @@ struct lanczos13m53
|
||||
static_cast<T>(0.1009141566987569892221439918230042368112e-8L),
|
||||
};
|
||||
T result = 0;
|
||||
T dz = z - 2;
|
||||
T z = dz + 2;
|
||||
for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
|
||||
{
|
||||
result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
|
||||
@@ -915,7 +915,7 @@ struct lanczos17m64
|
||||
}
|
||||
|
||||
template<class T>
|
||||
static T lanczos_sum_near_2(const T& z)
|
||||
static T lanczos_sum_near_2(const T& dz)
|
||||
{
|
||||
static const T d[16] = {
|
||||
static_cast<T>(23.56409085052261327114594781581930373708L),
|
||||
@@ -936,7 +936,7 @@ struct lanczos17m64
|
||||
static_cast<T>(-0.7398586479708476329563577384044188912075e-16L),
|
||||
};
|
||||
T result = 0;
|
||||
T dz = z - 2;
|
||||
T z = dz + 2;
|
||||
for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
|
||||
{
|
||||
result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
|
||||
@@ -1108,7 +1108,7 @@ struct lanczos24m113
|
||||
}
|
||||
|
||||
template<class T>
|
||||
static T lanczos_sum_near_2(const T& z)
|
||||
static T lanczos_sum_near_2(const T& dz)
|
||||
{
|
||||
static const T d[23] = {
|
||||
static_cast<T>(61.4165001061101455341808888883960361969557848005400286332291451422461117307237198559485365L),
|
||||
@@ -1136,7 +1136,7 @@ struct lanczos24m113
|
||||
static_cast<T>(0.410605371184590959139968810080063542546949719163227555918846829816144878123034347778284006e-25L),
|
||||
};
|
||||
T result = 0;
|
||||
T dz = z - 2;
|
||||
T z = dz + 2;
|
||||
for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
|
||||
{
|
||||
result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
|
||||
|
||||
@@ -3,7 +3,7 @@
|
||||
// Boost Software License, Version 1.0. (See accompanying file
|
||||
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#define L22
|
||||
//#define L22
|
||||
#include "../tools/ntl_rr_lanczos.hpp"
|
||||
#include "../tools/ntl_rr_digamma.hpp"
|
||||
#include <boost/math/tools/ntl.hpp>
|
||||
@@ -183,16 +183,6 @@ NTL::RR f(const NTL::RR& x, int variant)
|
||||
return l;
|
||||
}
|
||||
return expm1(x) / x;
|
||||
case 12:
|
||||
// log1p:
|
||||
if(x == 0)
|
||||
return 1;
|
||||
return log1p(x) / x;
|
||||
case 13:
|
||||
// log1p:
|
||||
if(x == 0)
|
||||
return 1;
|
||||
return log1p(-x) / (-x);
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
@@ -157,12 +157,16 @@ void test_spots(RealType T)
|
||||
static_cast<RealType>(9.999500016666250E-005L), // p
|
||||
static_cast<RealType>(1-9.999500016666250E-005L), //q
|
||||
tolerance);
|
||||
/*
|
||||
// This test data appears to be erroneous, MathCad appears
|
||||
// to suffer from cancellation error as x -> 0
|
||||
test_spot(
|
||||
static_cast<RealType>(1), // lambda
|
||||
static_cast<RealType>(0.0000001), // x
|
||||
static_cast<RealType>(9.999999499998730E-008L), // p
|
||||
static_cast<RealType>(1-9.999999499998730E-008L), //q
|
||||
tolerance);
|
||||
*/
|
||||
|
||||
BOOST_CHECK_CLOSE(
|
||||
::boost::math::pdf(
|
||||
@@ -285,18 +289,12 @@ Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_exponential_
|
||||
Running 1 test case...
|
||||
Tolerance for type float is 0.000596046 %
|
||||
Tolerance for type double is 1.11022e-012 %
|
||||
../../../../../../boost-sandbox/libs/math_functions/test/test_exponential_dist.cpp(43): error in "test_main_caller( argc, argv )": difference between ::boost::math::cdf( exponential_distribution<RealType>(l), x){9.9999995000000163e-008} and p{9.9999994999987298e-008} exceeds 1.11022e-012%
|
||||
../../../../../../boost-sandbox/libs/math_functions/test/test_exponential_dist.cpp(57): error in "test_main_caller( argc, argv )": difference between ::boost::math::quantile( exponential_distribution<RealType>(l), p){9.9999999999987131e-008} and x{9.9999999999999995e-008} exceeds 1.11022e-012%
|
||||
Tolerance for type long double is 1.11022e-012 %
|
||||
../../../../../../boost-sandbox/libs/math_functions/test/test_exponential_dist.cpp(43): error in "test_main_caller( argc, argv )": difference between ::boost::math::cdf( exponential_distribution<RealType>(l), x){9.9999995000000163e-008} and p{9.9999994999987298e-008} exceeds 1.11022e-012%
|
||||
../../../../../../boost-sandbox/libs/math_functions/test/test_exponential_dist.cpp(57): error in "test_main_caller( argc, argv )": difference between ::boost::math::quantile( exponential_distribution<RealType>(l), p){9.9999999999987131e-008} and x{9.9999999999999995e-008} exceeds 1.11022e-012%
|
||||
Tolerance for type class boost::math::concepts::real_concept is 1.11022e-012 %
|
||||
*** 4 failures detected in test suite "Test Program"
|
||||
Project : error PRJ0019: A tool returned an error code from "Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_exponential_dist.exe""
|
||||
Build Time 0:05
|
||||
Build log was saved at "file://i:\boost-06-05-03-1300\libs\math\test\Math_test\test_exponential\Debug\BuildLog.htm"
|
||||
test_exponential_dist - 5 error(s), 0 warning(s)
|
||||
========== Build: 0 succeeded, 1 failed, 0 up-to-date, 0 skipped ==========
|
||||
Build log was saved at "file://c:\data\boost\sandbox\boost-sandbox\libs\math_functions\IDE\math_toolkit\test_distribution\Debug\BuildLog.htm"
|
||||
test_distribution - 0 error(s), 1 warning(s)
|
||||
========== Build: 1 succeeded, 0 failed, 0 up-to-date, 0 skipped ==========
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -4199,7 +4199,7 @@ void print_code(const lanczos_info<T>& l, const char* name)
|
||||
|
||||
std::cout <<
|
||||
"\n template<class T>\n"
|
||||
" static T lanczos_sum_near_2(const T& z)\n"
|
||||
" static T lanczos_sum_near_2(const T& dz)\n"
|
||||
" {\n"
|
||||
" static const T d[" << l2.n-1 << "] = {\n";
|
||||
|
||||
@@ -4213,7 +4213,7 @@ void print_code(const lanczos_info<T>& l, const char* name)
|
||||
std::cout <<
|
||||
" };\n"
|
||||
" T result = 0;\n"
|
||||
" T dz = z - 2;\n"
|
||||
" T z = dz + 2;\n"
|
||||
" for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)\n"
|
||||
" {\n"
|
||||
" result += (-d[k-1]*dz)/(z + k*z + k*k - 1);\n"
|
||||
|
||||
@@ -111,7 +111,7 @@ struct lanczos13UDT
|
||||
}
|
||||
|
||||
template<class T>
|
||||
static T lanczos_sum_near_2(const T& z)
|
||||
static T lanczos_sum_near_2(const T& dz)
|
||||
{
|
||||
static const T d[12] = {
|
||||
boost::lexical_cast<T>("26.96979819614830698367887026728396466395"),
|
||||
@@ -128,7 +128,7 @@ struct lanczos13UDT
|
||||
boost::lexical_cast<T>("-0.9685385411006641478305219367315965391289e-9"),
|
||||
};
|
||||
T result = 0;
|
||||
T dz = z - 2;
|
||||
T z = dz + 2;
|
||||
for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
|
||||
{
|
||||
result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
|
||||
@@ -291,7 +291,7 @@ struct lanczos22UDT
|
||||
}
|
||||
|
||||
template<class T>
|
||||
static T lanczos_sum_near_2(const T& z)
|
||||
static T lanczos_sum_near_2(const T& dz)
|
||||
{
|
||||
static const T d[21] = {
|
||||
boost::lexical_cast<T>("75.39272007105208086018421070699575462226"),
|
||||
@@ -317,7 +317,7 @@ struct lanczos22UDT
|
||||
boost::lexical_cast<T>("0.6580207998808093935798753964580596673177e-19"),
|
||||
};
|
||||
T result = 0;
|
||||
T dz = z - 2;
|
||||
T z = dz + 2;
|
||||
for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
|
||||
{
|
||||
result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
|
||||
@@ -524,7 +524,7 @@ struct lanczos31UDT
|
||||
}
|
||||
|
||||
template<class T>
|
||||
static T lanczos_sum_near_2(const T& z)
|
||||
static T lanczos_sum_near_2(const T& dz)
|
||||
{
|
||||
static const T d[30] = {
|
||||
boost::lexical_cast<T>("147.9979641587472136175636384176549713358"),
|
||||
@@ -559,7 +559,7 @@ struct lanczos31UDT
|
||||
boost::lexical_cast<T>("-0.1060655106553177007425710511436497259484e-29"),
|
||||
};
|
||||
T result = 0;
|
||||
T dz = z - 2;
|
||||
T z = dz + 2;
|
||||
for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
|
||||
{
|
||||
result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
|
||||
@@ -802,7 +802,7 @@ struct lanczos61UDT
|
||||
}
|
||||
|
||||
template<class T>
|
||||
static T lanczos_sum_near_2(const T& z)
|
||||
static T lanczos_sum_near_2(const T& dz)
|
||||
{
|
||||
using namespace boost;
|
||||
static const T d[60] = {
|
||||
@@ -868,7 +868,7 @@ struct lanczos61UDT
|
||||
boost::lexical_cast<T>("-0.684593199208628857931267904308244537968349564351534581234005234847904343404822808648361291e-70"),
|
||||
};
|
||||
T result = 0;
|
||||
T dz = z - 2;
|
||||
T z = dz + 2;
|
||||
for(int k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
|
||||
{
|
||||
result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
|
||||
|
||||
Reference in New Issue
Block a user