Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve performance of the Complete Elliptic Integrals. #815

Merged
merged 7 commits into from
Aug 24, 2022
Merged

Improve performance of the Complete Elliptic Integrals. #815

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
d687773
Improve performance of the Complete Elliptic Integrals.
jzmaddock Aug 10, 2022
9302b79
Correct some accuracy issues in ellint specializations.
jzmaddock Aug 11, 2022
10bd7be
ellint_1/2 performance tweaks.
jzmaddock Aug 17, 2022
58ab7db
Define BOOST_FORCEINLINE in standalone mode.
jzmaddock Aug 18, 2022
7d219b2
Merge branch 'develop' into complete_elliptic_optimizations
jzmaddock Aug 22, 2022
c8fdd65
Re-factor nested calls from ellint_3->ellint_1.
jzmaddock Aug 23, 2022
de7928e
Disable new ellint code for clang-7.
jzmaddock Aug 24, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 22 additions & 22 deletions doc/graphs/elliptic_integral_e__80_bit_long_double.svg
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
73 changes: 21 additions & 52 deletions doc/graphs/elliptic_integral_e__double.svg
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
46 changes: 23 additions & 23 deletions doc/graphs/elliptic_integral_k__80_bit_long_double.svg
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
8 changes: 4 additions & 4 deletions doc/graphs/elliptic_integral_k__double.svg
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
18 changes: 9 additions & 9 deletions doc/graphs/plot_1d_errors.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ void plot_errors_1d(F f, Real start, Real end, unsigned points, const char* func
Real pos = start;
Real interval = (end - start) / points;

std::map<Real, Real> points_upper, points_lower;
std::map<double, double> points_upper, points_lower;

Real max_distance(0), min_distance(0), max_error(0), max_error_location(0);

Expand All @@ -60,19 +60,19 @@ void plot_errors_1d(F f, Real start, Real end, unsigned points, const char* func
Real exact_value = static_cast<Real>(f(mp_type(pos)));
Real distance = boost::math::sign(found_value - exact_value) * boost::math::epsilon_difference(found_value, exact_value);
Real bin = start + ((end - start) / num_bins) * boost::math::itrunc(num_bins * (pos - start) / (end - start));
if (points_lower.find(bin) == points_lower.end())
points_lower[bin] = 0;
if (points_upper.find(bin) == points_upper.end())
points_upper[bin] = 0;
if (points_lower.find((double)bin) == points_lower.end())
points_lower[(double)bin] = 0;
if (points_upper.find((double)bin) == points_upper.end())
points_upper[(double)bin] = 0;
if (distance > 0)
{
if (points_upper[bin] < distance)
points_upper[bin] = (std::min)(distance, max_y_scale);
if (points_upper[(double)bin] < distance)
points_upper[(double)bin] = (double)(std::min)(distance, max_y_scale);
}
else
{
if (points_lower[bin] > distance)
points_lower[bin] = (std::max)(distance, -max_y_scale);
if (points_lower[(double)bin] > distance)
points_lower[(double)bin] = (double)(std::max)(distance, -max_y_scale);
}
if (max_distance < distance)
max_distance = (std::min)(distance, max_y_scale);
Expand Down
14 changes: 12 additions & 2 deletions doc/sf/ellint_legendre.qbk
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,12 @@ NTL::RR at 1000-bit precision and this implementation.

[heading Implementation]

These functions are implemented in terms of Carlson's integrals using the relations:
For up to 80-bit long double precision the complete versions of these functions
are implemented as Taylor series expansions as in:
"Fast computation of complete elliptic integrals and Jacobian elliptic functions",
Celestial Mechanics and Dynamical Astronomy, April 2012.

Otherwise these functions are implemented in terms of Carlson's integrals using the relations:

[equation ellint19]

Expand Down Expand Up @@ -198,7 +203,12 @@ NTL::RR at 1000-bit precision and this implementation.

[heading Implementation]

These functions are implemented in terms of Carlson's integrals
For up to 80-bit long double precision the complete versions of these functions
are implemented as Taylor series expansions as in:
"Fast computation of complete elliptic integrals and Jacobian elliptic functions",
Celestial Mechanics and Dynamical Astronomy, April 2012.

Otherwise these functions are implemented in terms of Carlson's integrals
using the relations:

[equation ellint21]
Expand Down
26 changes: 13 additions & 13 deletions include/boost/math/policies/error_handling.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -751,7 +751,7 @@ namespace detail
{

template <class R, class T, class Policy>
inline bool check_overflow(T val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error))
BOOST_FORCEINLINE bool check_overflow(T val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error))
{
BOOST_MATH_STD_USING
if(fabs(val) > tools::max_value<R>())
Expand All @@ -763,7 +763,7 @@ inline bool check_overflow(T val, R* result, const char* function, const Policy&
return false;
}
template <class R, class T, class Policy>
inline bool check_overflow(std::complex<T> val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error))
BOOST_FORCEINLINE bool check_overflow(std::complex<T> val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error))
{
typedef typename R::value_type r_type;
r_type re, im;
Expand All @@ -773,7 +773,7 @@ inline bool check_overflow(std::complex<T> val, R* result, const char* function,
return r;
}
template <class R, class T, class Policy>
inline bool check_underflow(T val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error))
BOOST_FORCEINLINE bool check_underflow(T val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error))
{
if((val != 0) && (static_cast<R>(val) == 0))
{
Expand All @@ -783,7 +783,7 @@ inline bool check_underflow(T val, R* result, const char* function, const Policy
return false;
}
template <class R, class T, class Policy>
inline bool check_underflow(std::complex<T> val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error))
BOOST_FORCEINLINE bool check_underflow(std::complex<T> val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error))
{
typedef typename R::value_type r_type;
r_type re, im;
Expand All @@ -793,7 +793,7 @@ inline bool check_underflow(std::complex<T> val, R* result, const char* function
return r;
}
template <class R, class T, class Policy>
inline bool check_denorm(T val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error))
BOOST_FORCEINLINE bool check_denorm(T val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error))
{
BOOST_MATH_STD_USING
if((fabs(val) < static_cast<T>(tools::min_value<R>())) && (static_cast<R>(val) != 0))
Expand All @@ -804,7 +804,7 @@ inline bool check_denorm(T val, R* result, const char* function, const Policy& p
return false;
}
template <class R, class T, class Policy>
inline bool check_denorm(std::complex<T> val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error))
BOOST_FORCEINLINE bool check_denorm(std::complex<T> val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error))
{
typedef typename R::value_type r_type;
r_type re, im;
Expand All @@ -816,28 +816,28 @@ inline bool check_denorm(std::complex<T> val, R* result, const char* function, c

// Default instantiations with ignore_error policy.
template <class R, class T>
inline constexpr bool check_overflow(T /* val */, R* /* result */, const char* /* function */, const overflow_error<ignore_error>&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T))
BOOST_FORCEINLINE constexpr bool check_overflow(T /* val */, R* /* result */, const char* /* function */, const overflow_error<ignore_error>&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T))
{ return false; }
template <class R, class T>
inline constexpr bool check_overflow(std::complex<T> /* val */, R* /* result */, const char* /* function */, const overflow_error<ignore_error>&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T))
BOOST_FORCEINLINE constexpr bool check_overflow(std::complex<T> /* val */, R* /* result */, const char* /* function */, const overflow_error<ignore_error>&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T))
{ return false; }
template <class R, class T>
inline constexpr bool check_underflow(T /* val */, R* /* result */, const char* /* function */, const underflow_error<ignore_error>&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T))
BOOST_FORCEINLINE constexpr bool check_underflow(T /* val */, R* /* result */, const char* /* function */, const underflow_error<ignore_error>&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T))
{ return false; }
template <class R, class T>
inline constexpr bool check_underflow(std::complex<T> /* val */, R* /* result */, const char* /* function */, const underflow_error<ignore_error>&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T))
BOOST_FORCEINLINE constexpr bool check_underflow(std::complex<T> /* val */, R* /* result */, const char* /* function */, const underflow_error<ignore_error>&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T))
{ return false; }
template <class R, class T>
inline constexpr bool check_denorm(T /* val */, R* /* result*/, const char* /* function */, const denorm_error<ignore_error>&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T))
BOOST_FORCEINLINE constexpr bool check_denorm(T /* val */, R* /* result*/, const char* /* function */, const denorm_error<ignore_error>&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T))
{ return false; }
template <class R, class T>
inline constexpr bool check_denorm(std::complex<T> /* val */, R* /* result*/, const char* /* function */, const denorm_error<ignore_error>&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T))
BOOST_FORCEINLINE constexpr bool check_denorm(std::complex<T> /* val */, R* /* result*/, const char* /* function */, const denorm_error<ignore_error>&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T))
{ return false; }

} // namespace detail

template <class R, class Policy, class T>
inline R checked_narrowing_cast(T val, const char* function) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && is_noexcept_error_policy<Policy>::value)
BOOST_FORCEINLINE R checked_narrowing_cast(T val, const char* function) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && is_noexcept_error_policy<Policy>::value)
{
typedef typename Policy::overflow_error_type overflow_type;
typedef typename Policy::underflow_error_type underflow_type;
Expand Down
Loading