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

[Ellint Pi] Fix spurious underflow in sqrt(-v*N).

Browse Source
This commit is contained in:
jzmaddock
2014-12-27 13:30:51 +00:00
parent 47c2f9254c
commit e6b86ba93d
4 changed files with 24 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
include/boost/math/special_functions/ellint_3.hpp
View File

@@ -145,7 +145,11 @@ T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol)
T k2 = k * k;
T N = (k2 - v) / (1 - v);
T Nm1 = (1 - k2) / (1 - v);
T p2 = sqrt(-v * N);
T p2 = -v * N;
if(p2 <= tools::min_value<T>())
p2 = sqrt(-v) * sqrt(N);
else
p2 = sqrt(p2);
T delta = sqrt(1 - k2 * sphi * sphi);
if(N > k2)
{

7
test/float128/test_ellint_3.cpp
View File

@@ -22,6 +22,13 @@ void expected_results()
".*", // test type(s)
".*Large.*", // test data group
".*", 75, 40); // test function
add_expected_result(
".*", // compiler
".*", // stdlib
".*", // platform
".*", // test type(s)
".*Mathworld.*", // test data group
".*", 600, 200); // test function
add_expected_result(
".*", // compiler
".*", // stdlib

2
test/test_ellint_3.cpp
View File

@@ -69,7 +69,7 @@ void expected_results()
".*", // platform
largest_type, // test type(s)
".*Mathworld.*", // test data group
".*", 100, 50); // test function
".*", 600, 200); // test function
add_expected_result(
".*", // compiler
".*", // stdlib

17
test/test_ellint_3.hpp
View File

@@ -81,7 +81,7 @@ void test_spots(T, const char* type_name)
{
BOOST_MATH_STD_USING
// function values calculated on http://functions.wolfram.com/
static const boost::array<boost::array<T, 4>, 61> data1 = {{
static const boost::array<boost::array<T, 4>, 65> data1 = {{
{{ SC_(1.0), SC_(-1.0), SC_(0.0), SC_(-1.557407724654902230506974807458360173087) }},
{{ SC_(0.0), SC_(-4.0), SC_(0.4), SC_(-4.153623371196831087495427530365430979011) }},
{{ SC_(0.0), SC_(8.0), SC_(-0.6), SC_(8.935930619078575123490612395578518914416) }},
@@ -116,6 +116,10 @@ void test_spots(T, const char* type_name)
{ { std::numeric_limits<T>::max_exponent > 600 ? -ldexp(T(1), 604) : T(0), -ldexp(T(1), -816), -ldexp(T(1), -510), std::numeric_limits<T>::max_exponent > 600 ? SC_(-2.28835573409367516299079046268930870596307631872422530813192e-246) : SC_(-2.28835573409367516299079046268930870596307631872422530813192e-246) } },
{ { -ldexp(T(1), -622), -ldexp(T(1), -800), ldexp(T(1), -132), SC_(-1.49969681389563095481764443762806535353996169623910829793733e-241) } },
{ { -ldexp(T(1), -622), -ldexp(T(1), -800), -ldexp(T(1), -132), SC_(-1.49969681389563095481764443762806535353996169623910829793733e-241) } },
{ { -ldexp(T(1), -562), ldexp(T(1), -140), ldexp(T(1), -256), SC_(7.174648137343063403129495466444370592154941142407760751e-43) } },
{ { -ldexp(T(1), -562), -ldexp(T(1), -140), ldexp(T(1), -256), SC_(-7.17464813734306340312949546644437059215494114240776075e-43) } },
{ { ldexp(T(1), -688), -ldexp(T(1), -243), ldexp(T(1), -193), SC_(-7.07474928033336903711649944600608732865822749854620171e-74) } },
{ { -ldexp(T(1), -688), -ldexp(T(1), -243), ldexp(T(1), -193), SC_(-7.07474928033336903711649944600608732865822749854620171e-74) } },
// Special cases where k = 0:
{ { SC_(0.5), SC_(1.0), SC_(0.0), SC_(1.17881507892743738986863357869566288974084658835353613038547) } },
{ { SC_(-0.5), SC_(1.0), SC_(0.0), SC_(0.888286691263535380266337576823783210424994266596287990733270) } },
@@ -140,11 +144,12 @@ void test_spots(T, const char* type_name)
{ { SC_(1.0), SC_(-2.0), ldexp(T(1), -200), SC_(2.18503986326151899164330610231368254343201774622766316456295) } },
{ { SC_(1.0), ldexp(T(1.0), -150), ldexp(T(1), -200), SC_(7.006492321624085354618647916449580656401309709382578858e-46) } },
{ { SC_(1.0), -ldexp(T(1.0), -150), ldexp(T(1), -200), SC_(-7.006492321624085354618647916449580656401309709382578858e-46) } },
// Previosuly unsupported region with v > 1 and |phi| > PI/2:
{ { SC_(20.0), ldexp(T(1), -10) + (std::numeric_limits<T>::digits > 50 ? boost::math::constants::pi<T>() : 0), SC_(0.5), SC_(0.000976568747692583207373162769739923352941882247128764119469) } },
{ { SC_(20.0), -ldexp(T(1), -10) - (std::numeric_limits<T>::digits > 50 ? boost::math::constants::pi<T>() : 0), SC_(0.5), SC_(-0.000976568747692583207373162769739923352941882247128764119469) } },
{ { SC_(1.0) + ldexp(T(1), -6), (std::numeric_limits<T>::digits > 50 ? SC_(1.695796326794896619231321691639751442098584699687552910487472) : SC_(0.125)), SC_(0.5), (std::numeric_limits<T>::digits > 50 ? SC_(-27.1556829547879272134963658442519620188811279238098426583417) : SC_(0.125747519536592933239312389845380877095789502237158621392763)) } },
{ { SC_(1.0) + ldexp(T(1), -6), (std::numeric_limits<T>::digits > 50 ? SC_(-1.695796326794896619231321691639751442098584699687552910487472) : SC_(-0.125)), SC_(0.5), (std::numeric_limits<T>::digits > 50 ? SC_(27.1556829547879272134963658442519620188811279238098426583417) : SC_(-0.125747519536592933239312389845380877095789502237158621392763)) } },
// Previously unsupported region with v > 1 and |phi| > PI/2, this is the only region
// with high-ish error rates caused by argument reduction by Pi:
{ { SC_(20.0), ldexp(T(1647611), -19), SC_(0.5), SC_(0.000975940902692994840122139131147517258405256880370413541280) } },
{ { SC_(20.0), -ldexp(T(1647611), -19), SC_(0.5), SC_(-0.000975940902692994840122139131147517258405256880370413541280) } },
{ { SC_(1.0) + ldexp(T(1), -6), ldexp(T(889085), -19), SC_(0.5), SC_(-27.1647225624906589308619292363045712770651414487085887109197) } },
{ { SC_(1.0) + ldexp(T(1), -6), -ldexp(T(889085), -19), SC_(0.5), SC_(27.1647225624906589308619292363045712770651414487085887109197) } },
// Phi = 0:
{ { SC_(1.0), SC_(0.0), SC_(0.5), SC_(0.0) } },
{ { SC_(-1.0), SC_(0.0), SC_(0.5), SC_(0.0) } },