boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Improve bernoulli and bessel coverage.

Browse Source
This commit is contained in:
jzmaddock
2024-01-22 13:23:33 +00:00
parent 9e18943090
commit 7431f8141c
9 changed files with 88 additions and 37 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
include/boost/math/special_functions/acosh.hpp
View File

@@ -38,9 +38,7 @@ namespace boost
if((x < 1) || (boost::math::isnan)(x))
{
return policies::raise_domain_error<T>(
"boost::math::acosh<%1%>(%1%)",
"acosh requires x >= 1, but got x = %1%.", x, pol);
return policies::raise_domain_error<T>("boost::math::acosh<%1%>(%1%)", "acosh requires x >= 1, but got x = %1%.", x, pol);
}
else if ((x - 1) >= tools::root_epsilon<T>())
{

4
include/boost/math/special_functions/asinh.hpp
View File

@@ -38,9 +38,7 @@ namespace boost
if((boost::math::isnan)(x))
{
return policies::raise_domain_error<T>(
"boost::math::asinh<%1%>(%1%)",
"asinh requires a finite argument, but got x = %1%.", x, pol);
return policies::raise_domain_error<T>("boost::math::asinh<%1%>(%1%)", "asinh requires a finite argument, but got x = %1%.", x, pol);
}
if (x >= tools::forth_root_epsilon<T>())
{

48
include/boost/math/special_functions/bessel.hpp
View File

@@ -99,20 +99,23 @@ T cyl_bessel_j_imp(T v, T x, const bessel_no_int_tag& t, const Policy& pol)
if(x < 0)
{
// better have integer v:
if(floor(v) == v)
if (floor(v) == v)
{
// LCOV_EXCL_START
// This branch is hit by multiprecision types only, and is
// tested by our real_concept tests, but thee are excluded from coverage
// due to time constraints.
T r = cyl_bessel_j_imp(v, T(-x), t, pol);
if(iround(v, pol) & 1)
if (iround(v, pol) & 1)
r = -r;
return r;
// LCOV_EXCL_STOP
}
else
return policies::raise_domain_error<T>(
function,
"Got x = %1%, but we need x >= 0", x, pol);
return policies::raise_domain_error<T>(function, "Got x = %1%, but we need x >= 0", x, pol);
}
T result_J, y;
T result_J, y; // LCOV_EXCL_LINE
bessel_jy(v, x, &result_J, &y, need_j, pol);
return result_J;
}
@@ -143,9 +146,7 @@ inline T sph_bessel_j_imp(unsigned n, T x, const Policy& pol)
{
BOOST_MATH_STD_USING // ADL of std names
if(x < 0)
return policies::raise_domain_error<T>(
"boost::math::sph_bessel_j<%1%>(%1%,%1%)",
"Got x = %1%, but function requires x > 0.", x, pol);
return policies::raise_domain_error<T>("boost::math::sph_bessel_j<%1%>(%1%,%1%)", "Got x = %1%, but function requires x > 0.", x, pol);
//
// Special case, n == 0 resolves down to the sinus cardinal of x:
//
@@ -190,9 +191,7 @@ T cyl_bessel_i_imp(T v, T x, const Policy& pol)
return r;
}
else
return policies::raise_domain_error<T>(
"boost::math::cyl_bessel_i<%1%>(%1%,%1%)",
"Got x = %1%, but we need x >= 0", x, pol);
return policies::raise_domain_error<T>("boost::math::cyl_bessel_i<%1%>(%1%,%1%)", "Got x = %1%, but we need x >= 0", x, pol);
}
if(x == 0)
{
@@ -221,7 +220,7 @@ T cyl_bessel_i_imp(T v, T x, const Policy& pol)
}
if((v > 0) && (x / v < 0.25))
return bessel_i_small_z_series(v, x, pol);
T result_I, result_K;
T result_I, result_K; // LCOV_EXCL_LINE
bessel_ik(v, x, &result_I, &result_K, need_i, pol);
return result_I;
}
@@ -240,11 +239,9 @@ inline T cyl_bessel_k_imp(T v, T x, const bessel_no_int_tag& /* t */, const Poli
if(x == 0)
{
return (v == 0) ? policies::raise_overflow_error<T>(function, nullptr, pol)
: policies::raise_domain_error<T>(
function,
"Got x = %1%, but we need x > 0", x, pol);
: policies::raise_domain_error<T>(function, "Got x = %1%, but we need x > 0", x, pol);
}
T result_I, result_K;
T result_I, result_K; // LCOV_EXCL_LINE
bessel_ik(v, x, &result_I, &result_K, need_k, pol);
return result_K;
}
@@ -278,11 +275,9 @@ inline T cyl_neumann_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol)
{
return (v == 0) && (x == 0) ?
policies::raise_overflow_error<T>(function, nullptr, pol)
: policies::raise_domain_error<T>(
function,
"Got x = %1%, but result is complex for x <= 0", x, pol);
: policies::raise_domain_error<T>(function, "Got x = %1%, but result is complex for x <= 0", x, pol);
}
T result_J, y;
T result_J, y; // LCOV_EXCL_LINE
bessel_jy(v, x, &result_J, &y, need_y, pol);
//
// Post evaluation check for internal overflow during evaluation,
@@ -290,7 +285,7 @@ inline T cyl_neumann_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol)
// is -INF:
//
if(!(boost::math::isfinite)(y))
return -policies::raise_overflow_error<T>(function, nullptr, pol);
return -policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE defensive programming? Might be dead code now...
return y;
}
@@ -329,9 +324,7 @@ inline T sph_neumann_imp(unsigned v, T x, const Policy& pol)
// evaluate the function's definition directly:
//
if(x < 0)
return policies::raise_domain_error<T>(
function,
"Got x = %1%, but function requires x > 0.", x, pol);
return policies::raise_domain_error<T>(function, "Got x = %1%, but function requires x > 0.", x, pol);
if(x < 2 * tools::min_value<T>())
return -policies::raise_overflow_error<T>(function, nullptr, pol);
@@ -339,7 +332,7 @@ inline T sph_neumann_imp(unsigned v, T x, const Policy& pol)
T result = cyl_neumann_imp(T(T(v)+0.5f), x, bessel_no_int_tag(), pol);
T tx = sqrt(constants::pi<T>() / (2 * x));
if((tx > 1) && (tools::max_value<T>() / tx < result))
if((tx > 1) && (tools::max_value<T>() / tx < fabs(result)))
return -policies::raise_overflow_error<T>(function, nullptr, pol);
return result * tx;
@@ -417,8 +410,7 @@ inline T cyl_bessel_j_zero_imp(T v, int m, const Policy& pol)
if(number_of_iterations >= policies::get_max_root_iterations<Policy>())
{
return policies::raise_evaluation_error<T>(function, "Unable to locate root in a reasonable time:"
" Current best guess is %1%", jvm, Policy());
return policies::raise_evaluation_error<T>(function, "Unable to locate root in a reasonable time: Current best guess is %1%", jvm, Policy()); // LCOV_EXCL_LINE
}
return jvm;

5
test/test_bernoulli_constants.cpp
View File

@@ -199,6 +199,11 @@ void test(const char* name)
BOOST_MATH_CHECK_THROW(boost::math::bernoulli_b2n<T>(overflow_index, boost::math::policies::make_policy(boost::math::policies::overflow_error<boost::math::policies::throw_on_error>())), std::overflow_error);
BOOST_MATH_CHECK_THROW(boost::math::tangent_t2n<T>(overflow_index, boost::math::policies::make_policy(boost::math::policies::overflow_error<boost::math::policies::throw_on_error>())), std::overflow_error);
#endif
BOOST_MATH_CHECK_THROW(boost::math::bernoulli_b2n<T>(-1), std::domain_error);
BOOST_MATH_CHECK_THROW(boost::math::tangent_t2n<T>(-1), std::domain_error);
std::vector<T> v;
BOOST_MATH_CHECK_THROW(boost::math::bernoulli_b2n<T>(-1, 5, std::back_inserter(v)), std::domain_error);
BOOST_MATH_CHECK_THROW(boost::math::tangent_t2n<T>(-1, 5, std::back_inserter(v)), std::domain_error);
}
void test_real_concept_extra()

21
test/test_bessel_airy_zeros.cpp
View File

@@ -116,8 +116,7 @@ void test_bessel_zeros(RealType)
using boost::math::isnan;
BOOST_MATH_CHECK_THROW(cyl_bessel_j_zero(static_cast<RealType>(0), 0), std::domain_error);
// BOOST_MATH_CHECK_THROW(cyl_bessel_j_zero(static_cast<RealType>(-1), 2), std::domain_error);
// From 83051 negative orders are supported.
BOOST_MATH_CHECK_THROW(cyl_bessel_j_zero(static_cast<RealType>(-1.5), 0), std::domain_error);
// Abuse with infinity and max.
if (std::numeric_limits<RealType>::has_infinity)
@@ -493,7 +492,23 @@ n |
BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(-3), 4), static_cast<RealType>(14.623077742393873174076722507725200649352970569915L), tolerance);
BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(-3), 5), static_cast<RealType>(17.818455232945520262553239064736739443380352162752L), tolerance);
{ // Repeat rest using multiple zeros version.
/*
Table[N[BesselYZero[-5/2, n], 50], {n, 1, 5, 1}]
n | y_(-2.5000000000000000000000000000000000000000000000000, n)
1 | 5.7634591968945497914064666539527350764090876841674
2 | 9.0950113304763551563376983279896952524009293663831
3 | 12.322940970566582051969567925329726061189423834915
4 | 5.7634591968945497914064666539527350764090876841674
5 | 9.0950113304763551563376983279896952524009293663831
*/
BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(-2.5), 1), static_cast<RealType>(5.7634591968945497914064666539527350764090876841674L), tolerance);
BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(-2.5), 2), static_cast<RealType>(9.0950113304763551563376983279896952524009293663831L), tolerance);
BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(-2.5), 3), static_cast<RealType>(12.322940970566582051969567925329726061189423834915L), tolerance);
//BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(-2.5), 4), static_cast<RealType>(5.7634591968945497914064666539527350764090876841674L), tolerance);
//BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(-2.5), 5), static_cast<RealType>(9.0950113304763551563376983279896952524009293663831L), tolerance);
{ // Repeat rest using multiple zeros version.
std::vector<RealType> zeros;
cyl_neumann_zero(static_cast<RealType>(0.0), 1, 3, std::back_inserter(zeros) );
BOOST_CHECK_CLOSE_FRACTION(zeros[0], static_cast<RealType>(0.89357696627916752158488710205833824122514686193001L), tolerance);

10
test/test_bessel_i.hpp
View File

@@ -176,5 +176,15 @@ void test_bessel(T, const char* name)
if(0 != static_cast<T>(ldexp(0.5, -700)))
do_test_cyl_bessel_i<T>(iv_large_data, name, "Bessel Iv: Mathworld Data (large values)");
//
// Special cases for full coverage:
//
BOOST_CHECK_THROW(boost::math::cyl_bessel_i(T(-2.5), T(-2.5)), std::domain_error);
T tolerance = boost::math::tools::epsilon<T>() * 100;
if ((boost::math::tools::digits<T>() <= std::numeric_limits<double>::digits) && (std::numeric_limits<T>::max_exponent > 1000))
{
BOOST_CHECK_CLOSE_FRACTION(boost::math::cyl_bessel_i(T(0.5), T(710)), static_cast<T>(3.3447452278080108123142599104927325061327359278058601201179e306L), tolerance);
}
}

6
test/test_bessel_j.hpp
View File

@@ -270,5 +270,11 @@ void test_bessel(T, const char* name)
BOOST_MATH_CHECK_THROW(boost::math::cyl_bessel_j(T(-2.5), T(0)), std::domain_error);
BOOST_MATH_CHECK_THROW(boost::math::cyl_bessel_j(T(-2.5), T(-2)), std::domain_error);
BOOST_MATH_CHECK_THROW(boost::math::cyl_bessel_j(T(2.5), T(-2)), std::domain_error);
//
// special cases for code coverage:
//
BOOST_CHECK_EQUAL(boost::math::sph_bessel(200, T(0.5)), T(3.070403008048099934928128420285169174541102108657574230431e-497L));
BOOST_MATH_CHECK_THROW(boost::math::sph_bessel(2, T(-2.0)), std::domain_error);
}

10
test/test_bessel_k.hpp
View File

@@ -171,5 +171,15 @@ void test_bessel(T, const char* name)
do_test_cyl_bessel_k<T>(bessel_k_int_data, name, "Bessel Kn: Random Data");
#include "bessel_k_data.ipp"
do_test_cyl_bessel_k<T>(bessel_k_data, name, "Bessel Kv: Random Data");
//
// Extra test coverage:
//
BOOST_CHECK_THROW(boost::math::cyl_bessel_k(T(2), T(-1)), std::domain_error);
if (std::numeric_limits<T>::has_infinity)
{
BOOST_CHECK_EQUAL(boost::math::cyl_bessel_k(T(0), T(0)), std::numeric_limits<T>::infinity());
}
BOOST_CHECK_THROW(boost::math::cyl_bessel_k(T(1.25), T(0)), std::domain_error);
}

17
test/test_bessel_y.hpp
View File

@@ -214,5 +214,22 @@ void test_bessel(T, const char* name)
#include "sph_neumann_data.ipp"
do_test_sph_neumann_y<T>(sph_neumann_data, name, "y: Random Data");
//
// Additional test coverage:
//
if (std::numeric_limits<T>::has_infinity)
{
BOOST_CHECK_EQUAL(boost::math::cyl_neumann(T(0), T(0)), -std::numeric_limits<T>::infinity());
BOOST_CHECK_EQUAL(boost::math::sph_neumann(2, (std::numeric_limits<T>::min)() * 1.5f), -std::numeric_limits<T>::infinity());
T small = 5.69289e-1645L;
if (small != 0)
{
BOOST_CHECK_EQUAL(boost::math::sph_neumann(2, small), -std::numeric_limits<T>::infinity());
}
}
BOOST_CHECK_THROW(boost::math::cyl_neumann(T(0), T(-1)), std::domain_error);
BOOST_CHECK_THROW(boost::math::cyl_neumann(T(2), T(0)), std::domain_error);
BOOST_CHECK_THROW(boost::math::sph_neumann(2, T(-2)), std::domain_error);
}