Skip to content

Commit

Permalink
Merge pull request #815 from boostorg/complete_elliptic_optimizations
Browse files Browse the repository at this point in the history 
Improve performance of the Complete Elliptic Integrals.
  • Loading branch information
@jzmaddock
jzmaddock authored Aug 24, 2022
2 parents f111b6a + de7928e commit d2a9c9e
Show file tree
Hide file tree
Showing 14 changed files with 1,729 additions and 153 deletions.
44 changes: 22 additions & 22 deletions doc/graphs/elliptic_integral_e__80_bit_long_double.svg
View file
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
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
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
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
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
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
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

0 comments on commit d2a9c9e

Please sign in to comment.