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] Try and actually use the correct equation for negative v in Pi[v, k] !

Browse Source
This commit is contained in:
jzmaddock
2014-12-21 19:08:18 +00:00
parent b885aa757a
commit 2ab808f547
2 changed files with 7 additions and 30 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -259,39 +259,14 @@ T ellint_pi_imp(T v, T k, T vc, const Policy& pol)
if(v < 0)
{
// Apply A&S 17.7.17:
T k2 = k * k;
T N = (k2 - v) / (1 - v);
T Nm1 = (1 - k2) / (1 - v);
T p2 = sqrt(-v * (k2 - v) / (1 - v));
T result = 0;
if(k2 < N)
{
result = boost::math::detail::ellint_pi_imp(N, k, Nm1, pol);
result *= sqrt(Nm1 * (1 - k2 / N));
}
/* else result = 0;
Result so far must be zero if k^2/N == 1, note that
since k^2 <= 1 it is necessarily true that k^2/N <= 1
even though slight roundoff error may actually yield
values > 1
*/
if(k != 0)
{
if((k2 <= tools::min_value<T>()) || (p2 <= tools::min_value<T>()))
{
// We need to use logs to calculate k2 / p2:
T lk2 = 2 * log(fabs(k)) - (log(-v) + log(k2 - v) - log(1 - v)) / 2;
lk2 = exp(lk2);
if(lk2 != 0)
result += ellint_k_imp(k, pol) * lk2;
}
else
{
result += ellint_k_imp(k, pol) * k2 / p2;
}
}
result /= sqrt((1 - v) * (1 - k2 / v));
result = boost::math::detail::ellint_pi_imp(N, k, Nm1, pol);
result *= -v * (1 - k2) / ((1 - v) * (k2 - v));
result += ellint_k_imp(k, pol) * k2 / (k2 - v);
return result;
}

4
test/test_ellint_3.hpp
View File

@@ -125,7 +125,7 @@ void test_spots(T, const char* type_name)
do_test_ellint_pi3<T>(ellint_pi3_large_data, type_name, "Elliptic Integral PI: Large Random Data");
// function values calculated on http://functions.wolfram.com/
static const boost::array<boost::array<T, 3>, 14> data2 = {{
static const boost::array<boost::array<T, 3>, 16> data2 = {{
{{ SC_(0.0), SC_(0.2), SC_(1.586867847454166237308008033828114192951) }},
{{ SC_(0.0), SC_(0.4), SC_(1.639999865864511206865258329748601457626) }},
{{ SC_(0.0), SC_(0.0), SC_(1.57079632679489661923132169163975144209858469968755291048747) }},
@@ -141,6 +141,8 @@ void test_spots(T, const char* type_name)
{ { SC_(1e-170), SC_(-1E-164), SC_(1.57079632679489661923132169163975144209858469968755291048747) } },
{ { -1.5f * ldexp(T(1), -52), SC_(-0.9375), SC_(2.48840049140103464299631535211815755485846563527849342319632) } },
{ { -1.5f * ldexp(T(1), -52), SC_(0.9375), SC_(2.48840049140103464299631535211815755485846563527849342319632) } },
{ { ldexp(T(1), -560), ldexp(T(1), -165), SC_(1.57079632679489661923132169163975144209858469968756130722545) } },
{ { ldexp(T(1), -560), -ldexp(T(1), -165), SC_(1.57079632679489661923132169163975144209858469968754451374949) } },
} };
do_test_ellint_pi2<T>(data2, type_name, "Complete Elliptic Integral PI: Mathworld Data");