Boost-Geometry-Utils
view release on metacpan or search on metacpan
src/boost/graph/named_function_params.hpp view on Meta::CPAN
struct orig_to_copy_t { };
struct root_vertex_t { };
struct polling_t { };
struct lookahead_t { };
struct in_parallel_t { };
struct attractive_force_t { };
struct repulsive_force_t { };
struct force_pairs_t { };
struct cooling_t { };
struct vertex_displacement_t { };
struct iterations_t { };
struct diameter_range_t { };
struct learning_constant_range_t { };
struct vertices_equivalent_t { };
struct edges_equivalent_t { };
struct index_in_heap_map_t { };
struct max_priority_queue_t { };
#define BOOST_BGL_DECLARE_NAMED_PARAMS \
BOOST_BGL_ONE_PARAM_CREF(weight_map, edge_weight) \
BOOST_BGL_ONE_PARAM_CREF(weight_map2, edge_weight2) \
src/boost/graph/named_function_params.hpp view on Meta::CPAN
BOOST_BGL_ONE_PARAM_CREF(vertex_invariant2, vertex_invariant2) \
BOOST_BGL_ONE_PARAM_CREF(vertex_max_invariant, vertex_max_invariant) \
BOOST_BGL_ONE_PARAM_CREF(polling, polling) \
BOOST_BGL_ONE_PARAM_CREF(lookahead, lookahead) \
BOOST_BGL_ONE_PARAM_CREF(in_parallel, in_parallel) \
BOOST_BGL_ONE_PARAM_CREF(displacement_map, vertex_displacement) \
BOOST_BGL_ONE_PARAM_CREF(attractive_force, attractive_force) \
BOOST_BGL_ONE_PARAM_CREF(repulsive_force, repulsive_force) \
BOOST_BGL_ONE_PARAM_CREF(force_pairs, force_pairs) \
BOOST_BGL_ONE_PARAM_CREF(cooling, cooling) \
BOOST_BGL_ONE_PARAM_CREF(iterations, iterations) \
BOOST_BGL_ONE_PARAM_CREF(diameter_range, diameter_range) \
BOOST_BGL_ONE_PARAM_CREF(learning_constant_range, learning_constant_range) \
BOOST_BGL_ONE_PARAM_CREF(vertices_equivalent, vertices_equivalent) \
BOOST_BGL_ONE_PARAM_CREF(edges_equivalent, edges_equivalent) \
BOOST_BGL_ONE_PARAM_CREF(index_in_heap_map, index_in_heap_map) \
BOOST_BGL_ONE_PARAM_REF(max_priority_queue, max_priority_queue)
template <typename T, typename Tag, typename Base = no_property>
struct bgl_named_params
{
src/boost/math/policies/error_handling.hpp view on Meta::CPAN
return result;
if(detail::check_underflow<R>(val, &result, function, underflow_type()))
return result;
if(detail::check_denorm<R>(val, &result, function, denorm_type()))
return result;
return static_cast<R>(val);
}
template <class T, class Policy>
inline void check_series_iterations(const char* function, boost::uintmax_t max_iter, const Policy& pol)
{
if(max_iter >= policies::get_max_series_iterations<Policy>())
raise_evaluation_error<T>(
function,
"Series evaluation exceeded %1% iterations, giving up now.", static_cast<T>(static_cast<double>(max_iter)), pol);
}
template <class T, class Policy>
inline void check_root_iterations(const char* function, boost::uintmax_t max_iter, const Policy& pol)
{
if(max_iter >= policies::get_max_root_iterations<Policy>())
raise_evaluation_error<T>(
function,
"Root finding evaluation exceeded %1% iterations, giving up now.", static_cast<T>(static_cast<double>(max_iter)), pol);
}
} //namespace policies
#ifdef BOOST_MSVC
# pragma warning(pop)
#endif
}} // namespaces boost/math
src/boost/math/policies/policy.hpp view on Meta::CPAN
BOOST_MATH_META_INT(discrete_quantile_policy_type, discrete_quantile, BOOST_MATH_DISCRETE_QUANTILE_POLICY)
//
// Precision:
//
BOOST_MATH_META_INT(int, digits10, BOOST_MATH_DIGITS10_POLICY)
BOOST_MATH_META_INT(int, digits2, 0)
//
// Iterations:
//
BOOST_MATH_META_INT(unsigned long, max_series_iterations, BOOST_MATH_MAX_SERIES_ITERATION_POLICY)
BOOST_MATH_META_INT(unsigned long, max_root_iterations, BOOST_MATH_MAX_ROOT_ITERATION_POLICY)
//
// Define the names for each possible policy:
//
#define BOOST_MATH_PARAMETER(name)\
BOOST_PARAMETER_TEMPLATE_KEYWORD(name##_name)\
BOOST_PARAMETER_NAME(name##_name)
struct default_policy{};
namespace detail{
src/boost/math/policies/policy.hpp view on Meta::CPAN
typedef typename detail::find_arg<arg_list, is_promote_double<mpl::_1>, promote_double<> >::type promote_double_type;
//
// Discrete quantiles:
//
typedef typename detail::find_arg<arg_list, is_discrete_quantile<mpl::_1>, discrete_quantile<> >::type discrete_quantile_type;
//
// Mathematically undefined properties:
//
typedef typename detail::find_arg<arg_list, is_assert_undefined<mpl::_1>, assert_undefined<> >::type assert_undefined_type;
//
// Max iterations:
//
typedef typename detail::find_arg<arg_list, is_max_series_iterations<mpl::_1>, max_series_iterations<> >::type max_series_iterations_type;
typedef typename detail::find_arg<arg_list, is_max_root_iterations<mpl::_1>, max_root_iterations<> >::type max_root_iterations_type;
};
//
// These full specializations are defined to reduce the amount of
// template instantiations that have to take place when using the default
// policies, they have quite a large impact on compile times:
//
template <>
struct policy<default_policy, default_policy, default_policy, default_policy, default_policy, default_policy, default_policy, default_policy, default_policy, default_policy, default_policy>
{
public:
src/boost/math/policies/policy.hpp view on Meta::CPAN
typedef indeterminate_result_error<> indeterminate_result_error_type;
#if BOOST_MATH_DIGITS10_POLICY == 0
typedef digits2<> precision_type;
#else
typedef detail::precision<digits10<>, digits2<> >::type precision_type;
#endif
typedef promote_float<> promote_float_type;
typedef promote_double<> promote_double_type;
typedef discrete_quantile<> discrete_quantile_type;
typedef assert_undefined<> assert_undefined_type;
typedef max_series_iterations<> max_series_iterations_type;
typedef max_root_iterations<> max_root_iterations_type;
};
template <>
struct policy<detail::forwarding_arg1, detail::forwarding_arg2, default_policy, default_policy, default_policy, default_policy, default_policy, default_policy, default_policy, default_policy, default_policy>
{
public:
typedef domain_error<> domain_error_type;
typedef pole_error<> pole_error_type;
typedef overflow_error<> overflow_error_type;
typedef underflow_error<> underflow_error_type;
src/boost/math/policies/policy.hpp view on Meta::CPAN
typedef indeterminate_result_error<> indeterminate_result_error_type;
#if BOOST_MATH_DIGITS10_POLICY == 0
typedef digits2<> precision_type;
#else
typedef detail::precision<digits10<>, digits2<> >::type precision_type;
#endif
typedef promote_float<false> promote_float_type;
typedef promote_double<false> promote_double_type;
typedef discrete_quantile<> discrete_quantile_type;
typedef assert_undefined<> assert_undefined_type;
typedef max_series_iterations<> max_series_iterations_type;
typedef max_root_iterations<> max_root_iterations_type;
};
template <class Policy,
class A1 = default_policy,
class A2 = default_policy,
class A3 = default_policy,
class A4 = default_policy,
class A5 = default_policy,
class A6 = default_policy,
class A7 = default_policy,
src/boost/math/policies/policy.hpp view on Meta::CPAN
typedef typename detail::find_arg<arg_list, is_promote_double<mpl::_1>, typename Policy::promote_double_type >::type promote_double_type;
//
// Discrete quantiles:
//
typedef typename detail::find_arg<arg_list, is_discrete_quantile<mpl::_1>, typename Policy::discrete_quantile_type >::type discrete_quantile_type;
//
// Mathematically undefined properties:
//
typedef typename detail::find_arg<arg_list, is_assert_undefined<mpl::_1>, typename Policy::assert_undefined_type >::type assert_undefined_type;
//
// Max iterations:
//
typedef typename detail::find_arg<arg_list, is_max_series_iterations<mpl::_1>, typename Policy::max_series_iterations_type>::type max_series_iterations_type;
typedef typename detail::find_arg<arg_list, is_max_root_iterations<mpl::_1>, typename Policy::max_root_iterations_type>::type max_root_iterations_type;
//
// Define a typelist of the policies:
//
typedef mpl::vector<
domain_error_type,
pole_error_type,
overflow_error_type,
underflow_error_type,
denorm_error_type,
evaluation_error_type,
rounding_error_type,
indeterminate_result_error_type,
precision_type,
promote_float_type,
promote_double_type,
discrete_quantile_type,
assert_undefined_type,
max_series_iterations_type,
max_root_iterations_type> result_list;
//
// Remove all the policies that are the same as the default:
//
typedef typename mpl::remove_if<result_list, detail::is_default_policy<mpl::_> >::type reduced_list;
//
// Pad out the list with defaults:
//
typedef typename detail::append_N<reduced_list, default_policy, (14 - ::boost::mpl::size<reduced_list>::value)>::type result_type;
public:
typedef policy<
src/boost/math/policies/policy.hpp view on Meta::CPAN
} // namespace detail
template <class T, class Policy>
inline int digits(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T))
{
typedef mpl::bool_< std::numeric_limits<T>::is_specialized > tag_type;
return detail::digits_imp<T, Policy>(tag_type());
}
template <class Policy>
inline unsigned long get_max_series_iterations()
{
typedef typename Policy::max_series_iterations_type iter_type;
return iter_type::value;
}
template <class Policy>
inline unsigned long get_max_root_iterations()
{
typedef typename Policy::max_root_iterations_type iter_type;
return iter_type::value;
}
namespace detail{
template <class T, class Digits, class Small, class Default>
struct series_factor_calc
{
static T get()
{
src/boost/math/special_functions/bessel.hpp view on Meta::CPAN
unsigned v;
T mult;
T term;
};
template <class T, class Policy>
inline T sph_bessel_j_small_z_series(unsigned v, T x, const Policy& pol)
{
BOOST_MATH_STD_USING // ADL of std names
sph_bessel_j_small_z_series_term<T, Policy> s(v, x);
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
T zero = 0;
T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
#else
T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
#endif
policies::check_series_iterations<T>("boost::math::sph_bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
return result * sqrt(constants::pi<T>() / 4);
}
template <class T, class Policy>
T cyl_bessel_j_imp(T v, T x, const bessel_no_int_tag& t, const Policy& pol)
{
BOOST_MATH_STD_USING
static const char* function = "boost::math::bessel_j<%1%>(%1%,%1%)";
if(x < 0)
{
src/boost/math/special_functions/beta.hpp view on Meta::CPAN
}
}
else
{
// Non-normalised, just compute the power:
result = pow(x, a);
}
if(result < tools::min_value<T>())
return s0; // Safeguard: series can't cope with denorms.
ibeta_series_t<T> s(a, b, x, result);
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0);
policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (with lanczos)", max_iter, pol);
return result;
}
//
// Incomplete Beta series again, this time without Lanczos support:
//
template <class T, class Policy>
T ibeta_series(T a, T b, T x, T s0, const boost::math::lanczos::undefined_lanczos&, bool normalised, T* p_derivative, T y, const Policy& pol)
{
BOOST_MATH_STD_USING
src/boost/math/special_functions/beta.hpp view on Meta::CPAN
}
}
else
{
// Non-normalised, just compute the power:
result = pow(x, a);
}
if(result < tools::min_value<T>())
return s0; // Safeguard: series can't cope with denorms.
ibeta_series_t<T> s(a, b, x, result);
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0);
policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (without lanczos)", max_iter, pol);
return result;
}
//
// Continued fraction for the incomplete beta:
//
template <class T>
struct ibeta_fraction2_t
{
typedef std::pair<T, T> result_type;
src/boost/math/special_functions/beta.hpp view on Meta::CPAN
for(unsigned n = 1; n < sizeof(p)/sizeof(p[0]); ++n)
{
/*
// debugging code, enable this if you want to determine whether
// the table of Pn's is large enough...
//
static int max_count = 2;
if(n > max_count)
{
max_count = n;
std::cerr << "Max iterations in BGRAT was " << n << std::endl;
}
*/
//
// begin by evaluating the next Pn from Eq 9.4:
//
tnp1 += 2;
p[n] = 0;
T mbn = b - n;
unsigned tmp1 = 3;
for(unsigned m = 1; m < n; ++m)
src/boost/math/special_functions/detail/bessel_ik.hpp view on Meta::CPAN
}
else
{
prefix = v * log(x / 2) - boost::math::lgamma(v + 1, pol);
prefix = exp(prefix);
}
if(prefix == 0)
return prefix;
cyl_bessel_i_small_z<T, Policy> s(v, x);
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
T zero = 0;
T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
#else
T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
#endif
policies::check_series_iterations<T>("boost::math::bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
return prefix * result;
}
// Calculate K(v, x) and K(v+1, x) by method analogous to
// Temme, Journal of Computational Physics, vol 21, 343 (1976)
template <typename T, typename Policy>
int temme_ik(T v, T x, T* K, T* K1, const Policy& pol)
{
T f, h, p, q, coef, sum, sum1, tolerance;
T a, b, c, d, sigma, gamma1, gamma2;
src/boost/math/special_functions/detail/bessel_ik.hpp view on Meta::CPAN
BOOST_MATH_INSTRUMENT_VARIABLE(sigma);
BOOST_MATH_INSTRUMENT_CODE(sinh(sigma));
BOOST_MATH_INSTRUMENT_VARIABLE(gamma1);
BOOST_MATH_INSTRUMENT_VARIABLE(gamma2);
BOOST_MATH_INSTRUMENT_VARIABLE(c);
BOOST_MATH_INSTRUMENT_VARIABLE(d);
BOOST_MATH_INSTRUMENT_VARIABLE(a);
// series summation
tolerance = tools::epsilon<T>();
for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
{
f = (k * f + p + q) / (k*k - v*v);
p /= k - v;
q /= k + v;
h = p - k * f;
coef *= x * x / (4 * k);
sum += coef * f;
sum1 += coef * h;
if (abs(coef * f) < abs(sum) * tolerance)
{
break;
}
}
policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in temme_ik", k, pol);
*K = sum;
*K1 = 2 * sum1 / x;
return 0;
}
// Evaluate continued fraction fv = I_(v+1) / I_v, derived from
// Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73
template <typename T, typename Policy>
int CF1_ik(T v, T x, T* fv, const Policy& pol)
{
T C, D, f, a, b, delta, tiny, tolerance;
unsigned long k;
BOOST_MATH_STD_USING
// |x| <= |v|, CF1_ik converges rapidly
// |x| > |v|, CF1_ik needs O(|x|) iterations to converge
// modified Lentz's method, see
// Lentz, Applied Optics, vol 15, 668 (1976)
tolerance = 2 * tools::epsilon<T>();
BOOST_MATH_INSTRUMENT_VARIABLE(tolerance);
tiny = sqrt(tools::min_value<T>());
BOOST_MATH_INSTRUMENT_VARIABLE(tiny);
C = f = tiny; // b0 = 0, replace with tiny
D = 0;
for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
{
a = 1;
b = 2 * (v + k) / x;
C = b + a / C;
D = b + a * D;
if (C == 0) { C = tiny; }
if (D == 0) { D = tiny; }
D = 1 / D;
delta = C * D;
f *= delta;
BOOST_MATH_INSTRUMENT_VARIABLE(delta-1);
if (abs(delta - 1) <= tolerance)
{
break;
}
}
BOOST_MATH_INSTRUMENT_VARIABLE(k);
policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in CF1_ik", k, pol);
*fv = f;
return 0;
}
// Calculate K(v, x) and K(v+1, x) by evaluating continued fraction
// z1 / z0 = U(v+1.5, 2v+1, 2x) / U(v+0.5, 2v+1, 2x), see
// Thompson and Barnett, Computer Physics Communications, vol 47, 245 (1987)
template <typename T, typename Policy>
src/boost/math/special_functions/detail/bessel_ik.hpp view on Meta::CPAN
f = delta = D; // f1 = delta1 = D1, coincidence
prev = 0; // q0
current = 1; // q1
Q = C = -a; // Q1 = C1 because q1 = 1
S = 1 + Q * delta; // S1
BOOST_MATH_INSTRUMENT_VARIABLE(tolerance);
BOOST_MATH_INSTRUMENT_VARIABLE(a);
BOOST_MATH_INSTRUMENT_VARIABLE(b);
BOOST_MATH_INSTRUMENT_VARIABLE(D);
BOOST_MATH_INSTRUMENT_VARIABLE(f);
for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++) // starting from 2
{
// continued fraction f = z1 / z0
a -= 2 * (k - 1);
b += 2;
D = 1 / (b + a * D);
delta *= b * D - 1;
f += delta;
// series summation S = 1 + \sum_{n=1}^{\infty} C_n * z_n / z_0
q = (prev - (b - 2) * current) / a;
src/boost/math/special_functions/detail/bessel_ik.hpp view on Meta::CPAN
S += Q * delta;
// S converges slower than f
BOOST_MATH_INSTRUMENT_VARIABLE(Q * delta);
BOOST_MATH_INSTRUMENT_VARIABLE(abs(S) * tolerance);
if (abs(Q * delta) < abs(S) * tolerance)
{
break;
}
}
policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in CF2_ik", k, pol);
*Kv = sqrt(pi<T>() / (2 * x)) * exp(-x) / S;
*Kv1 = *Kv * (0.5f + v + x + (v * v - 0.25f) * f) / x;
BOOST_MATH_INSTRUMENT_VARIABLE(*Kv);
BOOST_MATH_INSTRUMENT_VARIABLE(*Kv1);
return 0;
}
enum{
src/boost/math/special_functions/detail/bessel_jy.hpp view on Meta::CPAN
h = p;
coef = 1;
sum = coef * g;
sum1 = coef * h;
T v2 = v * v;
T coef_mult = -x * x / 4;
// series summation
tolerance = policies::get_epsilon<T, Policy>();
for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
{
f = (k * f + p + q) / (k*k - v2);
p /= k - v;
q /= k + v;
g = f + e * q;
h = p - k * g;
coef *= coef_mult / k;
sum += coef * g;
sum1 += coef * h;
if (abs(coef * g) < abs(sum) * tolerance)
{
break;
}
}
policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in temme_jy", k, pol);
*Y = -sum;
*Y1 = -2 * sum1 / x;
return 0;
}
// Evaluate continued fraction fv = J_(v+1) / J_v, see
// Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73
template <typename T, typename Policy>
int CF1_jy(T v, T x, T* fv, int* sign, const Policy& pol)
{
T C, D, f, a, b, delta, tiny, tolerance;
unsigned long k;
int s = 1;
BOOST_MATH_STD_USING
// |x| <= |v|, CF1_jy converges rapidly
// |x| > |v|, CF1_jy needs O(|x|) iterations to converge
// modified Lentz's method, see
// Lentz, Applied Optics, vol 15, 668 (1976)
tolerance = 2 * policies::get_epsilon<T, Policy>();;
tiny = sqrt(tools::min_value<T>());
C = f = tiny; // b0 = 0, replace with tiny
D = 0;
for (k = 1; k < policies::get_max_series_iterations<Policy>() * 100; k++)
{
a = -1;
b = 2 * (v + k) / x;
C = b + a / C;
D = b + a * D;
if (C == 0) { C = tiny; }
if (D == 0) { D = tiny; }
D = 1 / D;
delta = C * D;
f *= delta;
if (D < 0) { s = -s; }
if (abs(delta - 1) < tolerance)
{ break; }
}
policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF1_jy", k / 100, pol);
*fv = -f;
*sign = s; // sign of denominator
return 0;
}
//
// This algorithm was originally written by Xiaogang Zhang
// using std::complex to perform the complex arithmetic.
// However, that turns out to 10x or more slower than using
// all real-valued arithmetic, so it's been rewritten using
src/boost/math/special_functions/detail/bessel_jy.hpp view on Meta::CPAN
if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
temp = Dr * Dr + Di * Di;
Dr = Dr / temp;
Di = -Di / temp;
delta_r = Cr * Dr - Ci * Di;
delta_i = Ci * Dr + Cr * Di;
temp = fr;
fr = temp * delta_r - fi * delta_i;
fi = temp * delta_i + fi * delta_r;
for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++)
{
a = k - 0.5f;
a *= a;
a -= v2;
bi += 2;
temp = Cr * Cr + Ci * Ci;
Cr = br + a * Cr / temp;
Ci = bi - a * Ci / temp;
Dr = br + a * Dr;
Di = bi + a * Di;
src/boost/math/special_functions/detail/bessel_jy.hpp view on Meta::CPAN
Dr = Dr / temp;
Di = -Di / temp;
delta_r = Cr * Dr - Ci * Di;
delta_i = Ci * Dr + Cr * Di;
temp = fr;
fr = temp * delta_r - fi * delta_i;
fi = temp * delta_i + fi * delta_r;
if (fabs(delta_r - 1) + fabs(delta_i) < tolerance)
break;
}
policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol);
*p = fr;
*q = fi;
return 0;
}
static const int need_j = 1;
static const int need_y = 2;
// Compute J(v, x) and Y(v, x) simultaneously by Steed's method, see
src/boost/math/special_functions/detail/bessel_jy_asym.hpp view on Meta::CPAN
T q = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, v) / boost::math::tgamma(1 + v, pol);
T f = (p - q) / v;
T g_prefix = boost::math::sin_pi(v / 2, pol);
g_prefix *= g_prefix * 2 / v;
T g = f + g_prefix * q;
T h = p;
T c_mult = -x * x / 4;
T y(c * g), y1(c * h);
for(int k = 1; k < policies::get_max_series_iterations<Policy>(); ++k)
{
f = (k * f + p + q) / (k*k - v*v);
p /= k - v;
q /= k + v;
c *= c_mult / k;
T c1 = pow(-x * x / 4, k) / factorial<T>(k, pol);
g = f + g_prefix * q;
h = -k * g + p;
y += c * g;
y1 += c * h;
src/boost/math/special_functions/detail/bessel_jy_series.hpp view on Meta::CPAN
}
else
{
prefix = v * log(x / 2) - boost::math::lgamma(v+1, pol);
prefix = exp(prefix);
}
if(0 == prefix)
return prefix;
bessel_j_small_z_series_term<T, Policy> s(v, x);
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
T zero = 0;
T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
#else
T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
#endif
policies::check_series_iterations<T>("boost::math::bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
return prefix * result;
}
template <class T, class Policy>
struct bessel_y_small_z_series_term_a
{
typedef T result_type;
bessel_y_small_z_series_term_a(T v_, T x)
: N(0), v(v_)
src/boost/math/special_functions/detail/bessel_jy_series.hpp view on Meta::CPAN
prefix -= log(tools::max_value<T>() / 4);
scale /= (tools::max_value<T>() / 4);
if(tools::log_max_value<T>() < prefix)
{
return -policies::raise_overflow_error<T>(function, 0, pol);
}
}
prefix = -exp(prefix);
}
bessel_y_small_z_series_term_a<T, Policy> s(v, x);
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
*pscale = scale;
#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
T zero = 0;
T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
#else
T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
#endif
policies::check_series_iterations<T>("boost::math::bessel_y_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
result *= prefix;
if(!need_logs)
{
prefix = boost::math::tgamma(-v, pol) * boost::math::cos_pi(v) * p / constants::pi<T>();
}
else
{
int s;
prefix = boost::math::lgamma(-v, &s, pol) + p;
prefix = exp(prefix) * s / constants::pi<T>();
}
bessel_y_small_z_series_term_b<T, Policy> s2(v, x);
max_iter = policies::get_max_series_iterations<Policy>();
#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
#else
T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
#endif
result -= scale * prefix * b;
return result;
}
template <class T, class Policy>
src/boost/math/special_functions/detail/erf_inv.hpp view on Meta::CPAN
// Generic version, get a guess that's accurate to 64-bits (10^-19)
//
T guess = erf_inv_imp(p, q, pol, static_cast<mpl::int_<64> const*>(0));
T result;
//
// If T has more bit's than 64 in it's mantissa then we need to iterate,
// otherwise we can just return the result:
//
if(policies::digits<T, Policy>() > 64)
{
boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
if(p <= 0.5)
{
result = tools::halley_iterate(detail::erf_roots<typename remove_cv<T>::type, Policy>(p, 1), guess, static_cast<T>(0), tools::max_value<T>(), (policies::digits<T, Policy>() * 2) / 3, max_iter);
}
else
{
result = tools::halley_iterate(detail::erf_roots<typename remove_cv<T>::type, Policy>(q, -1), guess, static_cast<T>(0), tools::max_value<T>(), (policies::digits<T, Policy>() * 2) / 3, max_iter);
}
policies::check_root_iterations<T>("boost::math::erf_inv<%1%>", max_iter, pol);
}
else
{
result = guess;
}
return result;
}
template <class T, class Policy>
struct erf_inv_initializer
src/boost/math/special_functions/detail/gamma_inva.hpp view on Meta::CPAN
}
else if(z > 0.5)
{
guess = z * 1.2f;
}
else
{
guess = -0.4f / log(z);
}
//
// Max iterations permitted:
//
boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
//
// Use our generic derivative-free root finding procedure.
// We could use Newton steps here, taking the PDF of the
// Poisson distribution as our derivative, but that's
// even worse performance-wise than the generic method :-(
//
std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, false, tol, max_iter, pol);
if(max_iter >= policies::get_max_root_iterations<Policy>())
policies::raise_evaluation_error<T>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol);
return (r.first + r.second) / 2;
}
} // namespace detail
template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
gamma_p_inva(T1 x, T2 p, const Policy& pol)
{
typedef typename tools::promote_args<T1, T2>::type result_type;
src/boost/math/special_functions/detail/ibeta_inv_ab.hpp view on Meta::CPAN
}
else
{
guess = (std::min)(T(b / 2), T(10));
}
}
else
factor = (v < sqrt(tools::epsilon<T>())) ? 2 : (guess < 20 ? 1.2f : 1.1f);
BOOST_MATH_INSTRUMENT_CODE("guess = " << guess);
//
// Max iterations permitted:
//
boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, swap_ab ? true : false, tol, max_iter, pol);
if(max_iter >= policies::get_max_root_iterations<Policy>())
policies::raise_evaluation_error<T>("boost::math::ibeta_invab_imp<%1%>(%1%,%1%,%1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol);
return (r.first + r.second) / 2;
}
} // namespace detail
template <class RT1, class RT2, class RT3, class Policy>
typename tools::promote_args<RT1, RT2, RT3>::type
ibeta_inva(RT1 b, RT2 x, RT3 p, const Policy& pol)
{
typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
src/boost/math/special_functions/detail/ibeta_inverse.hpp view on Meta::CPAN
// that the distribution is *almost* symetric so 1-x will be in the right
// ball park for the solution:
//
if((x - s_2) * eta < 0)
x = 1 - x;
#ifdef BOOST_INSTRUMENT
std::cout << "Estimating x with Temme method 2 (large eta): " << x << std::endl;
#endif
}
//
// The final step is a few Newton-Raphson iterations to
// clean up our approximation for x, this is pretty cheap
// in general, and very cheap compared to an incomplete beta
// evaluation. The limits set on x come from the observation
// that the sign of eta and x-sin^2(theta) are the same.
//
T lower, upper;
if(eta < 0)
{
lower = 0;
upper = s_2;
src/boost/math/special_functions/detail/ibeta_inverse.hpp view on Meta::CPAN
// may be very skewed, these are not related by x ~ 1-x we used when
// implementing section 3 above. However we know that:
//
// cross < x <= 1 ; iff eta < mu
// x == cross ; iff eta == mu
// 0 <= x < cross ; iff eta > mu
//
// Where cross == 1 / (1 + mu)
// Many thanks to Prof Temme for clarifying this point.
//
// Therefore we'll just jump straight into Newton iterations
// to solve Eq 4.2 using these bounds, and simple bisection
// as the first guess, in practice this converges pretty quickly
// and we only need a few digits correct anyway:
//
if(eta <= 0)
eta = tools::min_value<T>();
T u = eta - mu * log(eta) + (1 + mu) * log(1 + mu) - mu;
T cross = 1 / (1 + mu);
T lower = eta < mu ? cross : 0;
T upper = eta < mu ? 1 : cross;
src/boost/math/special_functions/detail/ibeta_inverse.hpp view on Meta::CPAN
// term to estimate x, but when we expect x to be large then
// this greatly underestimates x and leaves us trying to
// iterate "round the corner" which may take almost forever...
//
// We could use Temme's inverse gamma function case in that case,
// this works really rather well (albeit expensively) even though
// strictly speaking we're outside it's defined range.
//
// However it's expensive to compute, and an alternative approach
// which models the curve as a distorted quarter circle is much
// cheaper to compute, and still keeps the number of iterations
// required down to a reasonable level. With thanks to Prof Temme
// for this suggestion.
//
if(b < a)
{
std::swap(a, b);
std::swap(p, q);
invert = !invert;
}
if(pow(p, 1/a) < 0.5)
src/boost/math/special_functions/detail/ibeta_inverse.hpp view on Meta::CPAN
invert = !invert;
T l = 1 - upper;
T u = 1 - lower;
lower = l;
upper = u;
}
//
// lower bound for our search:
//
// We're not interested in denormalised answers as these tend to
// these tend to take up lots of iterations, given that we can't get
// accurate derivatives in this area (they tend to be infinite).
//
if(lower == 0)
{
if(invert && (py == 0))
{
//
// We're not interested in answers smaller than machine epsilon:
//
lower = boost::math::tools::epsilon<T>();
src/boost/math/special_functions/detail/ibeta_inverse.hpp view on Meta::CPAN
// Figure out how many digits to iterate towards:
//
int digits = boost::math::policies::digits<T, Policy>() / 2;
if((x < 1e-50) && ((a < 1) || (b < 1)))
{
//
// If we're in a region where the first derivative is very
// large, then we have to take care that the root-finder
// doesn't terminate prematurely. We'll bump the precision
// up to avoid this, but we have to take care not to set the
// precision too high or the last few iterations will just
// thrash around and convergence may be slow in this case.
// Try 3/4 of machine epsilon:
//
digits *= 3;
digits /= 2;
}
//
// Now iterate, we can use either p or q as the target here
// depending on which is smaller:
//
boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
x = boost::math::tools::halley_iterate(
boost::math::detail::ibeta_roots<T, Policy>(a, b, (p < q ? p : q), (p < q ? false : true)), x, lower, upper, digits, max_iter);
policies::check_root_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%)", max_iter, pol);
//
// We don't really want these asserts here, but they are useful for sanity
// checking that we have the limits right, uncomment if you suspect bugs *only*.
//
//BOOST_ASSERT(x != upper);
//BOOST_ASSERT((x != lower) || (x == boost::math::tools::min_value<T>()) || (x == boost::math::tools::epsilon<T>()));
//
// Tidy up, if we "lower" was too high then zero is the best answer we have:
//
if(x == lower)
src/boost/math/special_functions/detail/igamma_inverse.hpp view on Meta::CPAN
else
{
digits /= 2;
digits -= 1;
}
if((a < 0.125) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
digits = policies::digits<T, Policy>() - 2;
//
// Go ahead and iterate:
//
boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
guess = tools::halley_iterate(
detail::gamma_p_inverse_func<T, Policy>(a, p, false),
guess,
lower,
tools::max_value<T>(),
digits,
max_iter);
policies::check_root_iterations<T>(function, max_iter, pol);
BOOST_MATH_INSTRUMENT_VARIABLE(guess);
if(guess == lower)
guess = policies::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
return guess;
}
template <class T, class Policy>
T gamma_q_inv_imp(T a, T q, const Policy& pol)
{
BOOST_MATH_STD_USING // ADL of std functions.
src/boost/math/special_functions/detail/igamma_inverse.hpp view on Meta::CPAN
else
{
digits /= 2;
digits -= 1;
}
if((a < 0.125) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
digits = policies::digits<T, Policy>();
//
// Go ahead and iterate:
//
boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
guess = tools::halley_iterate(
detail::gamma_p_inverse_func<T, Policy>(a, q, true),
guess,
lower,
tools::max_value<T>(),
digits,
max_iter);
policies::check_root_iterations<T>(function, max_iter, pol);
if(guess == lower)
guess = policies::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
return guess;
}
} // namespace detail
template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
gamma_p_inv(T1 a, T2 p, const Policy& pol)
src/boost/math/special_functions/ellint_rc.hpp view on Meta::CPAN
if (2 * abs(S) < tolerance)
break;
T sx = sqrt(x);
T sy = sqrt(y);
lambda = 2 * sx * sy + y;
x = (x + lambda) / 4;
y = (y + lambda) / 4;
++k;
}while(k < policies::get_max_series_iterations<Policy>());
// Check to see if we gave up too soon:
policies::check_series_iterations<T>(function, k, pol);
// Taylor series expansion to the 5th order
value = (1 + S * S * (T(3) / 10 + S * (T(1) / 7 + S * (T(3) / 8 + S * T(9) / 22)))) / sqrt(u);
return value * prefix;
}
} // namespace detail
template <class T1, class T2, class Policy>
src/boost/math/special_functions/ellint_rd.hpp view on Meta::CPAN
T sy = sqrt(y);
T sz = sqrt(z);
lambda = sy * (sx + sz) + sz * sx; //sqrt(x * y) + sqrt(y * z) + sqrt(z * x);
sigma += factor / (sz * (z + lambda));
factor /= 4;
x = (x + lambda) / 4;
y = (y + lambda) / 4;
z = (z + lambda) / 4;
++k;
}
while(k < policies::get_max_series_iterations<Policy>());
// Check to see if we gave up too soon:
policies::check_series_iterations<T>(function, k, pol);
// Taylor series expansion to the 5th order
EA = X * Y;
EB = Z * Z;
EC = EA - EB;
ED = EA - 6 * EB;
EE = ED + EC + EC;
S1 = ED * (ED * T(9) / 88 - Z * EE * T(9) / 52 - T(3) / 14);
S2 = Z * (EE / 6 + Z * (-EC * T(9) / 22 + Z * EA * T(3) / 26));
value = 3 * sigma + factor * (1 + S1 + S2) / (u * sqrt(u));
src/boost/math/special_functions/ellint_rf.hpp view on Meta::CPAN
T sx = sqrt(x);
T sy = sqrt(y);
T sz = sqrt(z);
lambda = sy * (sx + sz) + sz * sx;
x = (x + lambda) / 4;
y = (y + lambda) / 4;
z = (z + lambda) / 4;
++k;
}
while(k < policies::get_max_series_iterations<Policy>());
// Check to see if we gave up too soon:
policies::check_series_iterations<T>(function, k, pol);
BOOST_MATH_INSTRUMENT_VARIABLE(k);
// Taylor series expansion to the 5th order
E2 = X * Y - Z * Z;
E3 = X * Y * Z;
value = (1 + E2*(E2/24 - E3*T(3)/44 - T(0.1)) + E3/14) / sqrt(u);
BOOST_MATH_INSTRUMENT_VARIABLE(value);
return value;
}
src/boost/math/special_functions/ellint_rj.hpp view on Meta::CPAN
alpha *= alpha;
beta = p * (p + lambda) * (p + lambda);
sigma += factor * boost::math::ellint_rc(alpha, beta, pol);
factor /= 4;
x = (x + lambda) / 4;
y = (y + lambda) / 4;
z = (z + lambda) / 4;
p = (p + lambda) / 4;
++k;
}
while(k < policies::get_max_series_iterations<Policy>());
// Check to see if we gave up too soon:
policies::check_series_iterations<T>(function, k, pol);
// Taylor series expansion to the 5th order
EA = X * Y + Y * Z + Z * X;
EB = X * Y * Z;
EC = P * P;
E2 = EA - 3 * EC;
E3 = EB + 2 * P * (EA - EC);
S1 = 1 + E2 * (E2 * T(9) / 88 - E3 * T(9) / 52 - T(3) / 14);
S2 = EB * (T(1) / 6 + P * (T(-6) / 22 + P * T(3) / 26));
S3 = P * ((EA - EC) / 3 - P * EA * T(3) / 22);
src/boost/math/special_functions/erf.hpp view on Meta::CPAN
return -erf_imp(T(-z), invert, pol, t);
else
return 1 + erf_imp(T(-z), false, pol, t);
}
T result;
if(!invert && (z > detail::erf_asymptotic_limit<T, Policy>()))
{
detail::erf_asympt_series_t<T> s(z);
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
result = boost::math::tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, 1);
policies::check_series_iterations<T>("boost::math::erf<%1%>(%1%, %1%)", max_iter, pol);
}
else
{
T x = z * z;
if(x < 0.6)
{
// Compute P:
result = z * exp(-x);
result /= sqrt(boost::math::constants::pi<T>());
if(result != 0)
src/boost/math/special_functions/expint.hpp view on Meta::CPAN
T b;
int i;
unsigned n;
};
template <class T, class Policy>
inline T expint_as_fraction(unsigned n, T z, const Policy& pol)
{
BOOST_MATH_STD_USING
BOOST_MATH_INSTRUMENT_VARIABLE(z)
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
expint_fraction<T> f(n, z);
T result = tools::continued_fraction_b(
f,
boost::math::policies::get_epsilon<T, Policy>(),
max_iter);
policies::check_series_iterations<T>("boost::math::expint_continued_fraction<%1%>(unsigned,%1%)", max_iter, pol);
BOOST_MATH_INSTRUMENT_VARIABLE(result)
BOOST_MATH_INSTRUMENT_VARIABLE(max_iter)
result = exp(-z) / result;
BOOST_MATH_INSTRUMENT_VARIABLE(result)
return result;
}
template <class T>
struct expint_series
{
src/boost/math/special_functions/expint.hpp view on Meta::CPAN
T z;
T x_k;
T denom;
T fact;
};
template <class T, class Policy>
inline T expint_as_series(unsigned n, T z, const Policy& pol)
{
BOOST_MATH_STD_USING
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
BOOST_MATH_INSTRUMENT_VARIABLE(z)
T result = 0;
T x_k = -1;
T denom = T(1) - n;
T fact = 1;
unsigned k = 0;
for(; k < n - 1;)
{
src/boost/math/special_functions/expint.hpp view on Meta::CPAN
x_k *= -z;
fact *= ++k;
}
BOOST_MATH_INSTRUMENT_VARIABLE(result)
result += pow(-z, static_cast<T>(n - 1))
* (boost::math::digamma(static_cast<T>(n)) - log(z)) / fact;
BOOST_MATH_INSTRUMENT_VARIABLE(result)
expint_series<T> s(k, z, x_k, denom, fact);
result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, result);
policies::check_series_iterations<T>("boost::math::expint_series<%1%>(unsigned,%1%)", max_iter, pol);
BOOST_MATH_INSTRUMENT_VARIABLE(result)
BOOST_MATH_INSTRUMENT_VARIABLE(max_iter)
return result;
}
template <class T, class Policy, class Tag>
T expint_imp(unsigned n, T z, const Policy& pol, const Tag& tag)
{
BOOST_MATH_STD_USING
static const char* function = "boost::math::expint<%1%>(unsigned, %1%)";
src/boost/math/special_functions/expint.hpp view on Meta::CPAN
T z;
};
template <class T, class Policy>
T expint_i_as_series(T z, const Policy& pol)
{
BOOST_MATH_STD_USING
T result = log(z); // (log(z) - log(1 / z)) / 2;
result += constants::euler<T>();
expint_i_series<T> s(z);
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, result);
policies::check_series_iterations<T>("boost::math::expint_i_series<%1%>(%1%)", max_iter, pol);
return result;
}
template <class T, class Policy, class Tag>
T expint_i_imp(T z, const Policy& pol, const Tag& tag)
{
static const char* function = "boost::math::expint<%1%>(%1%)";
if(z < 0)
return -expint_imp(1, T(-z), pol, tag);
if(z == 0)
src/boost/math/special_functions/expm1.hpp view on Meta::CPAN
{
if(x > 0)
return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
return -1;
}
return exp(x) - T(1);
}
if(a < tools::epsilon<T>())
return x;
detail::expm1_series<T> s(x);
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) && !BOOST_WORKAROUND(__EDG_VERSION__, <= 245)
T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter);
#else
T zero = 0;
T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, zero);
#endif
policies::check_series_iterations<T>("boost::math::expm1<%1%>(%1%)", max_iter, pol);
return result;
}
template <class T, class P>
T expm1_imp(T x, const mpl::int_<53>&, const P& pol)
{
BOOST_MATH_STD_USING
T a = fabs(x);
if(a > T(0.5L))
src/boost/math/special_functions/gamma.hpp view on Meta::CPAN
}
};
template <class T, class Policy>
inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0)
{
// Multiply result by ((z^a) * (e^-z) / a) to get the full
// lower incomplete integral. Then divide by tgamma(a)
// to get the normalised value.
lower_incomplete_gamma_series<T> s(a, z);
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
T factor = policies::get_epsilon<T, Policy>();
T result = boost::math::tools::sum_series(s, factor, max_iter, init_value);
policies::check_series_iterations<T>("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
return result;
}
//
// Fully generic tgamma and lgamma use the incomplete partial
// sums added together:
//
template <class T, class Policy>
T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l)
{
src/boost/math/special_functions/gamma.hpp view on Meta::CPAN
// Compute the full upper fraction (Q) when a is very small:
//
T result;
result = boost::math::tgamma1pm1(a, pol);
if(pgam)
*pgam = (result + 1) / a;
T p = boost::math::powm1(x, a, pol);
result -= p;
result /= a;
detail::small_gamma2_series<T> s(a, x);
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>() - 10;
p += 1;
if(pderivative)
*pderivative = p / (*pgam * exp(x));
T init_value = invert ? *pgam : 0;
result = -p * tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, (init_value - result) / p);
policies::check_series_iterations<T>("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
if(invert)
result = -result;
return result;
}
//
// Upper gamma fraction for integer a:
//
template <class T, class Policy>
inline T finite_gamma_q(T a, T x, Policy const& pol, T* pderivative = 0)
{
src/boost/math/special_functions/log1p.hpp view on Meta::CPAN
function, 0, pol);
result_type a = abs(result_type(x));
if(a > result_type(0.5f))
return log(1 + result_type(x));
// Note that without numeric_limits specialisation support,
// epsilon just returns zero, and our "optimisation" will always fail:
if(a < tools::epsilon<result_type>())
return x;
detail::log1p_series<result_type> s(x);
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) && !BOOST_WORKAROUND(__EDG_VERSION__, <= 245)
result_type result = tools::sum_series(s, policies::get_epsilon<result_type, Policy>(), max_iter);
#else
result_type zero = 0;
result_type result = tools::sum_series(s, policies::get_epsilon<result_type, Policy>(), max_iter, zero);
#endif
policies::check_series_iterations<T>(function, max_iter, pol);
return result;
}
template <class T, class Policy>
T log1p_imp(T const& x, const Policy& pol, const mpl::int_<53>&)
{ // The function returns the natural logarithm of 1 + x.
BOOST_MATH_STD_USING
static const char* function = "boost::math::log1p<%1%>(%1%)";
src/boost/math/special_functions/log1p.hpp view on Meta::CPAN
result_type a = abs(result_type(x));
if(a > result_type(0.95f))
return log(1 + result_type(x)) - result_type(x);
// Note that without numeric_limits specialisation support,
// epsilon just returns zero, and our "optimisation" will always fail:
if(a < tools::epsilon<result_type>())
return -x * x / 2;
boost::math::detail::log1p_series<T> s(x);
s();
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
T zero = 0;
T result = boost::math::tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, zero);
#else
T result = boost::math::tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter);
#endif
policies::check_series_iterations<T>(function, max_iter, pol);
return result;
}
template <class T>
inline typename tools::promote_args<T>::type log1pmx(T x)
{
return log1pmx(x, policies::policy<>());
}
} // namespace math
src/boost/math/special_functions/owens_t.hpp view on Meta::CPAN
while( true )
{
RealType term = ai*yi;
val += term;
if((yi != 0) && (fabs(val * lim) > fabs(term)))
break;
ii += 2;
yi = (1.0-hs*yi) / static_cast<RealType>(ii);
ai *= as;
if(ii > (std::min)(1500, (int)policies::get_max_series_iterations<Policy>()))
policies::raise_evaluation_error("boost::math::owens_t<%1%>", 0, val, pol);
} // while( true )
return val;
} // arg_type owens_t_T4(const arg_type h, const arg_type a, const unsigned short m)
// This routine dispatches the call to one of six subroutines, depending on the values
// of h and a.
// preconditions: h >= 0, 0<=a<=1, ah=a*h
src/boost/math/special_functions/owens_t.hpp view on Meta::CPAN
try
{
have_t2 = true;
p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
if(p2.second < target_precision)
return p2.first;
}
catch(const boost::math::evaluation_error&){} // T2 may fail and throw, that's OK
}
//
// If we haven't tried T1 yet, do it now - sometimes it succeeds and the number of iterations
// is fairly low compared to T4.
//
if(!have_t1)
{
try
{
have_t1 = true;
p1 = owens_t_T1_accelerated(h, a, forwarding_policy());
if(p1.second < target_precision)
return p1.first;
}
catch(const boost::math::evaluation_error&){} // T1 may fail and throw, that's OK
}
//
// If we haven't tried T2 yet, do it now - sometimes it succeeds and the number of iterations
// is fairly low compared to T4.
//
if(!have_t2)
{
try
{
have_t2 = true;
p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
if(p2.second < target_precision)
return p2.first;
src/boost/math/special_functions/zeta.hpp view on Meta::CPAN
return pow(static_cast<T>(k++), s);
}
private:
T s;
unsigned k;
};
template <class T, class Policy>
inline T zeta_series2_imp(T s, const Policy& pol)
{
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();;
zeta_series2<T> f(s);
T result = tools::sum_series(
f,
policies::get_epsilon<T, Policy>(),
max_iter);
policies::check_series_iterations<T>("boost::math::zeta_series2<%1%>(%1%)", max_iter, pol);
return result;
}
#endif
template <class T, class Policy>
T zeta_polynomial_series(T s, T sc, Policy const &)
{
//
// This is algorithm 3 from:
//
src/boost/math/special_functions/zeta.hpp view on Meta::CPAN
T zeta_imp_prec(T s, T sc, const Policy& pol, const mpl::int_<0>&)
{
BOOST_MATH_STD_USING
T result;
result = zeta_polynomial_series(s, sc, pol);
#if 0
// Old code archived for future reference:
//
// Only use power series if it will converge in 100
// iterations or less: the more iterations it consumes
// the slower convergence becomes so we have to be very
// careful in it's usage.
//
if (s > -log(tools::epsilon<T>()) / 4.5)
result = detail::zeta_series2_imp(s, pol);
else
result = detail::zeta_series_imp(s, sc, pol);
#endif
return result;
}
src/boost/math/tools/roots.hpp view on Meta::CPAN
max_iter -= count;
#ifdef BOOST_MATH_INSTRUMENT
std::cout << "Bisection iteration, final count = " << max_iter << std::endl;
static boost::uintmax_t max_count = 0;
if(max_iter > max_count)
{
max_count = max_iter;
std::cout << "Maximum iterations: " << max_iter << std::endl;
}
#endif
return std::make_pair(min, max);
}
template <class F, class T, class Tol>
inline std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter)
{
return bisect(f, min, max, tol, max_iter, policies::policy<>());
src/boost/math/tools/roots.hpp view on Meta::CPAN
max_iter -= count;
#ifdef BOOST_MATH_INSTRUMENT
std::cout << "Newton Raphson iteration, final count = " << max_iter << std::endl;
static boost::uintmax_t max_count = 0;
if(max_iter > max_count)
{
max_count = max_iter;
std::cout << "Maximum iterations: " << max_iter << std::endl;
}
#endif
return result;
}
template <class F, class T>
inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits)
{
boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
src/boost/math/tools/roots.hpp view on Meta::CPAN
max_iter -= count;
#ifdef BOOST_MATH_INSTRUMENT
std::cout << "Schroeder iteration, final count = " << max_iter << std::endl;
static boost::uintmax_t max_count = 0;
if(max_iter > max_count)
{
max_count = max_iter;
std::cout << "Maximum iterations: " << max_iter << std::endl;
}
#endif
return result;
}
template <class F, class T>
inline T schroeder_iterate(F f, T guess, T min, T max, int digits)
{
boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
src/boost/math/tools/toms748_solve.hpp view on Meta::CPAN
{
//
// Zero is to the right of b, so walk upwards
// until we find it:
//
while((boost::math::sign)(fb) == (boost::math::sign)(fa))
{
if(count == 0)
policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol);
//
// Heuristic: every 20 iterations we double the growth factor in case the
// initial guess was *really* bad !
//
if((max_iter - count) % 20 == 0)
factor *= 2;
//
// Now go ahead and move our guess by "factor":
//
a = b;
fa = fb;
b *= factor;
src/boost/math/tools/toms748_solve.hpp view on Meta::CPAN
if(fabs(a) < tools::min_value<T>())
{
// Escape route just in case the answer is zero!
max_iter -= count;
max_iter += 1;
return a > 0 ? std::make_pair(T(0), T(a)) : std::make_pair(T(a), T(0));
}
if(count == 0)
policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol);
//
// Heuristic: every 20 iterations we double the growth factor in case the
// initial guess was *really* bad !
//
if((max_iter - count) % 20 == 0)
factor *= 2;
//
// Now go ahead and move are guess by "factor":
//
b = a;
fb = fa;
a /= factor;
( run in 1.052 second using v1.01-cache-2.11-cpan-96521ef73a4 )