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_1/2 performance tweaks.

Browse Source 
Add Google bench to probe changes more easily.
Update graphs and docs.
This commit is contained in:
jzmaddock
2022-08-17 19:21:03 +01:00
parent 9302b79221
commit 10bd7be0f1
10 changed files with 566 additions and 198 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
doc/graphs/elliptic_integral_e__80_bit_long_double.svg
View File

File diff suppressed because one or more lines are too long
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 21 KiB

After

Width:  |  Height:  |  Size: 21 KiB

73
doc/graphs/elliptic_integral_e__double.svg
View File

File diff suppressed because one or more lines are too long
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 28 KiB

After

Width:  |  Height:  |  Size: 21 KiB

46
doc/graphs/elliptic_integral_k__80_bit_long_double.svg
View File

File diff suppressed because one or more lines are too long
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 21 KiB

After

Width:  |  Height:  |  Size: 21 KiB

8
doc/graphs/elliptic_integral_k__double.svg
View File

File diff suppressed because one or more lines are too long
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 21 KiB

After

Width:  |  Height:  |  Size: 21 KiB

18
doc/graphs/plot_1d_errors.cpp
View File

@@ -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);
@@ -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);

14
doc/sf/ellint_legendre.qbk
View File

@@ -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]
@@ -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]

26
include/boost/math/policies/error_handling.hpp
View File

@@ -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>())
@@ -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;
@@ -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))
{
@@ -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;
@@ -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))
@@ -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;
@@ -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;

62
include/boost/math/special_functions/ellint_1.hpp
View File

@@ -187,23 +187,11 @@ T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 2> const&)
// archived in the code below), but was found to have slightly higher error rates.
//
template <typename T, typename Policy>
T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 0> const&)
BOOST_FORCEINLINE T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 0> const&)
{
using std::abs;
using namespace boost::math::tools;
static const char* function = "boost::math::ellint_k<%1%>(%1%)";
if (abs(k) > 1)
{
return policies::raise_domain_error<T>(function,
"Got k = %1%, function requires |k| <= 1", k, pol);
}
if (abs(k) == 1)
{
return policies::raise_overflow_error<T>(function, nullptr, pol);
}
T m = k * k;
switch (static_cast<int>(m * 20))
@@ -448,6 +436,10 @@ T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 0> const&)
return boost::math::tools::evaluate_polynomial(coef, m - 0.875);
}
default:
//
// This handles all cases where m > 0.9,
// including all error handling:
//
return ellint_k_imp(k, pol, std::integral_constant<int, 2>());
#if 0
else
@@ -468,23 +460,11 @@ T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 0> const&)
}
}
template <typename T, typename Policy>
T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 1> const&)
BOOST_FORCEINLINE T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 1> const&)
{
using std::abs;
using namespace boost::math::tools;
static const char* function = "boost::math::ellint_k<%1%>(%1%)";
if (abs(k) > 1)
{
return policies::raise_domain_error<T>(function,
"Got k = %1%, function requires |k| <= 1", k, pol);
}
if (abs(k) == 1)
{
return policies::raise_overflow_error<T>(function, nullptr, pol);
}
T m = k * k;
switch (static_cast<int>(m * 20))
{
@@ -757,28 +737,16 @@ T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 1> const&)
return boost::math::tools::evaluate_polynomial(coef, m - 0.875L);
}
default:
//
// All cases where m > 0.9
// including all error handling:
//
return ellint_k_imp(k, pol, std::integral_constant<int, 2>());
#if 0
default:
{
T lambda_prime = (1 - sqrt(k)) / (2 * (1 + sqrt(k)));
T k_prime = ellint_k(sqrt((1 - k) * (1 + k))); // K(m')
T lambda_prime_4th = boost::math::pow<4>(lambda_prime);
T q_prime = ((((((20910 * lambda_prime_4th) + 1707) * lambda_prime_4th + 150) * lambda_prime_4th + 15) * lambda_prime_4th + 2) * lambda_prime_4th + 1) * lambda_prime;
/*T q_prime_2 = lambda_prime
+ 2 * boost::math::pow<5>(lambda_prime)
+ 15 * boost::math::pow<9>(lambda_prime)
+ 150 * boost::math::pow<13>(lambda_prime)
+ 1707 * boost::math::pow<17>(lambda_prime)
+ 20910 * boost::math::pow<21>(lambda_prime);*/
return -log(q_prime) * k_prime / boost::math::constants::pi<T>();
}
#endif
}
}
template <typename T, typename Policy>
inline typename tools::promote_args<T>::type ellint_1(T k, const Policy& pol, const std::true_type&)
BOOST_FORCEINLINE typename tools::promote_args<T>::type ellint_1(T k, const Policy& pol, const std::true_type&)
{
typedef typename tools::promote_args<T>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
@@ -790,7 +758,7 @@ inline typename tools::promote_args<T>::type ellint_1(T k, const Policy& pol, co
}
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const std::false_type&)
BOOST_FORCEINLINE typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const std::false_type&)
{
return boost::math::ellint_1(k, phi, policies::policy<>());
}
@@ -799,14 +767,14 @@ inline typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const s
// Complete elliptic integral (Legendre form) of the first kind
template <typename T>
inline typename tools::promote_args<T>::type ellint_1(T k)
BOOST_FORCEINLINE typename tools::promote_args<T>::type ellint_1(T k)
{
return ellint_1(k, policies::policy<>());
}
// Elliptic integral (Legendre form) of the first kind
template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const Policy& pol)
BOOST_FORCEINLINE typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const Policy& pol)
{
typedef typename tools::promote_args<T1, T2>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
@@ -814,7 +782,7 @@ inline typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const P
}
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi)
BOOST_FORCEINLINE typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi)
{
typedef typename policies::is_policy<T2>::type tag_type;
return detail::ellint_1(k, phi, tag_type());

41
include/boost/math/special_functions/ellint_2.hpp
View File

@@ -189,21 +189,11 @@ T ellint_e_imp(T k, const Policy& pol, std::integral_constant<int, 2> const&)
// existing routines.
//
template <typename T, typename Policy>
T ellint_e_imp(T k, const Policy& pol, std::integral_constant<int, 0> const&)
BOOST_FORCEINLINE T ellint_e_imp(T k, const Policy& pol, std::integral_constant<int, 0> const&)
{
using std::abs;
using namespace boost::math::tools;
if (abs(k) > 1)
{
return policies::raise_domain_error<T>("boost::math::ellint_e<%1%>(%1%)",
"Got k = %1%, function requires |k| <= 1", k, pol);
}
if (abs(k) == 1)
{
return static_cast<T>(1);
}
T m = k * k;
switch (static_cast<int>(20 * m))
{
@@ -430,24 +420,19 @@ T ellint_e_imp(T k, const Policy& pol, std::integral_constant<int, 0> const&)
return boost::math::tools::evaluate_polynomial(coef, m - 0.875);
}
default:
//
// All cases where m > 0.9
// including all error handling:
//
return ellint_e_imp(k, pol, std::integral_constant<int, 2>());
}
}
template <typename T, typename Policy>
T ellint_e_imp(T k, const Policy& pol, std::integral_constant<int, 1> const&)
BOOST_FORCEINLINE T ellint_e_imp(T k, const Policy& pol, std::integral_constant<int, 1> const&)
{
using std::abs;
using namespace boost::math::tools;
if (abs(k) > 1)
{
return policies::raise_domain_error<T>("boost::math::ellint_e<%1%>(%1%)",
"Got k = %1%, function requires |k| <= 1", k, pol);
}
if (abs(k) == 1)
{
return static_cast<T>(1);
}
T m = k * k;
switch (static_cast<int>(20 * m))
{
@@ -708,12 +693,16 @@ T ellint_e_imp(T k, const Policy& pol, std::integral_constant<int, 1> const&)
return boost::math::tools::evaluate_polynomial(coef, m - 0.875L);
}
default:
//
// All cases where m > 0.9
// including all error handling:
//
return ellint_e_imp(k, pol, std::integral_constant<int, 2>());
}
}
template <typename T, typename Policy>
inline typename tools::promote_args<T>::type ellint_2(T k, const Policy& pol, const std::true_type&)
BOOST_FORCEINLINE typename tools::promote_args<T>::type ellint_2(T k, const Policy& pol, const std::true_type&)
{
typedef typename tools::promote_args<T>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
@@ -726,7 +715,7 @@ inline typename tools::promote_args<T>::type ellint_2(T k, const Policy& pol, co
// Elliptic integral (Legendre form) of the second kind
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const std::false_type&)
BOOST_FORCEINLINE typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const std::false_type&)
{
return boost::math::ellint_2(k, phi, policies::policy<>());
}
@@ -735,21 +724,21 @@ inline typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const s
// Complete elliptic integral (Legendre form) of the second kind
template <typename T>
inline typename tools::promote_args<T>::type ellint_2(T k)
BOOST_FORCEINLINE typename tools::promote_args<T>::type ellint_2(T k)
{
return ellint_2(k, policies::policy<>());
}
// Elliptic integral (Legendre form) of the second kind
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi)
BOOST_FORCEINLINE typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi)
{
typedef typename policies::is_policy<T2>::type tag_type;
return detail::ellint_2(k, phi, tag_type());
}
template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const Policy& pol)
BOOST_FORCEINLINE typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const Policy& pol)
{
typedef typename tools::promote_args<T1, T2>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;

432
reporting/performance/ellint_performance.cpp Normal file
View File

@@ -0,0 +1,432 @@
// (C) Copyright John Maddock 2022.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
//
// The purpose of this performance test is to probe the performance
// the of polynomial approximation to Elliptic integrals K and E.
//
// In the original paper:
// "Fast computation of complete elliptic integrals and Jacobian elliptic functions",
// Celestial Mechanics and Dynamical Astronomy, April 2012.
// they claim performance comparable to std::log, so we add comparison to log here
// as a reference "known good".
//
// We also measure the effect of disabling overflow errors, and of a "bare-bones"
// implementation which lacks a lot of our usual boilerplate.
//
// Note in the performance test routines, we had to sum the results of the function calls
// otherwise some msvc results were unrealistically fast, despite the use of
// benchmark::DoNotOptimize.
//
// Some sample results from msvc, and taking log(double) as a score of 1:
//
// ellint_1_performance<double> 1.7
// ellint_1_performance_no_overflow_check<double> 1.6
// ellint_2_performance<double> 1.8
// ellint_2_performance_no_overflow_check<double> 1.6
// basic_ellint_rational_performance<double> 1.6
//
// We can in fact get basic_ellint_rational_performance to much the same performance as log
// ONLY if we remove all error hangling for cases with m > 0.9. In particular the code appears
// to be ultra-sensitive to the presence of "if" statements which significantly hamper optimisation.
//
// Performance with gcc-cygwin appears to be broasly similar.
//
#include <vector>
#include <benchmark/benchmark.h>
#include <boost/math/special_functions/ellint_1.hpp>
#include <boost/math/special_functions/ellint_2.hpp>
#include <boost/math/tools/random_vector.hpp>
using boost::math::generate_random_uniform_vector;
template <typename Real>
const std::vector<Real>& get_test_data()
{
static const std::vector<Real> data = generate_random_uniform_vector<Real>(1024, 12345, 0.0, 0.9);
return data;
}
template <typename Real>
Real& get_result_data()
{
static Real data = 0;
return data;
}
template <class T>
BOOST_FORCEINLINE T basic_ellint_rational(T k)
{
using namespace boost::math::tools;
static const char* function = "boost::math::ellint_k<%1%>(%1%)";
T m = k * k;
switch (static_cast<int>(m * 20))
{
case 0:
case 1:
//if (m < 0.1)
{
constexpr T coef[] =
{
1.591003453790792180,
0.416000743991786912,
0.245791514264103415,
0.179481482914906162,
0.144556057087555150,
0.123200993312427711,
0.108938811574293531,
0.098853409871592910,
0.091439629201749751,
0.085842591595413900,
0.081541118718303215,
0.078199656811256481910
};
return boost::math::tools::evaluate_polynomial(coef, m - 0.05);
}
case 2:
case 3:
//else if (m < 0.2)
{
constexpr T coef[] =
{
1.635256732264579992,
0.471190626148732291,
0.309728410831499587,
0.252208311773135699,
0.226725623219684650,
0.215774446729585976,
0.213108771877348910,
0.216029124605188282,
0.223255831633057896,
0.234180501294209925,
0.248557682972264071,
0.266363809892617521
};
return boost::math::tools::evaluate_polynomial(coef, m - 0.15);
}
case 4:
case 5:
//else if (m < 0.3)
{
constexpr T coef[] =
{
1.685750354812596043,
0.541731848613280329,
0.401524438390690257,
0.369642473420889090,
0.376060715354583645,
0.405235887085125919,
0.453294381753999079,
0.520518947651184205,
0.609426039204995055,
0.724263522282908870,
0.871013847709812357,
1.057652872753547036
};
return boost::math::tools::evaluate_polynomial(coef, m - 0.25);
}
case 6:
case 7:
//else if (m < 0.4)
{
constexpr T coef[] =
{
1.744350597225613243,
0.634864275371935304,
0.539842564164445538,
0.571892705193787391,
0.670295136265406100,
0.832586590010977199,
1.073857448247933265,
1.422091460675497751,
1.920387183402304829,
2.632552548331654201,
3.652109747319039160,
5.115867135558865806,
7.224080007363877411
};
return boost::math::tools::evaluate_polynomial(coef, m - 0.35);
}
case 8:
case 9:
//else if (m < 0.5)
{
constexpr T coef[] =
{
1.813883936816982644,
0.763163245700557246,
0.761928605321595831,
0.951074653668427927,
1.315180671703161215,
1.928560693477410941,
2.937509342531378755,
4.594894405442878062,
7.330071221881720772,
11.87151259742530180,
19.45851374822937738,
32.20638657246426863,
53.73749198700554656,
90.27388602940998849
};
return boost::math::tools::evaluate_polynomial(coef, m - 0.45);
}
case 10:
case 11:
//else if (m < 0.6)
{
constexpr T coef[] =
{
1.898924910271553526,
0.950521794618244435,
1.151077589959015808,
1.750239106986300540,
2.952676812636875180,
5.285800396121450889,
9.832485716659979747,
18.78714868327559562,
36.61468615273698145,
72.45292395127771801,
145.1079577347069102,
293.4786396308497026,
598.3851815055010179,
1228.420013075863451,
2536.529755382764488
};
return boost::math::tools::evaluate_polynomial(coef, m - 0.55);
}
case 12:
case 13:
//else if (m < 0.7)
{
constexpr T coef[] =
{
2.007598398424376302,
1.248457231212347337,
1.926234657076479729,
3.751289640087587680,
8.119944554932045802,
18.66572130873555361,
44.60392484291437063,
109.5092054309498377,
274.2779548232413480,
697.5598008606326163,
1795.716014500247129,
4668.381716790389910,
12235.76246813664335,
32290.17809718320818,
85713.07608195964685,
228672.1890493117096,
612757.2711915852774
};
return boost::math::tools::evaluate_polynomial(coef, m - 0.65);
}
case 14:
case 15:
//else if (m < 0.8)
{
constexpr T coef[] =
{
2.156515647499643235,
1.791805641849463243,
3.826751287465713147,
10.38672468363797208,
31.40331405468070290,
100.9237039498695416,
337.3268282632272897,
1158.707930567827917,
4060.990742193632092,
14454.00184034344795,
52076.66107599404803,
189493.6591462156887,
695184.5762413896145,
2567994.048255284686,
9541921.966748386322,
35634927.44218076174,
133669298.4612040871,
503352186.6866284541,
1901975729.538660119,
7208915015.330103756
};
return boost::math::tools::evaluate_polynomial(coef, m - 0.75);
}
case 16:
//else if (m < 0.85)
{
constexpr T coef[] =
{
2.318122621712510589,
2.616920150291232841,
7.897935075731355823,
30.50239715446672327,
131.4869365523528456,
602.9847637356491617,
2877.024617809972641,
14110.51991915180325,
70621.44088156540229,
358977.2665825309926,
1847238.263723971684,
9600515.416049214109,
50307677.08502366879,
265444188.6527127967,
1408862325.028702687,
7515687935.373774627
};
return boost::math::tools::evaluate_polynomial(coef, m - 0.825);
}
case 17:
//else if (m < 0.90)
{
constexpr T coef[] =
{
2.473596173751343912,
3.727624244118099310,
15.60739303554930496,
84.12850842805887747,
506.9818197040613935,
3252.277058145123644,
21713.24241957434256,
149037.0451890932766,
1043999.331089990839,
7427974.817042038995,
53503839.67558661151,
389249886.9948708474,
2855288351.100810619,
21090077038.76684053,
156699833947.7902014,
1170222242422.439893,
8777948323668.937971,
66101242752484.95041,
499488053713388.7989,
37859743397240299.20
};
return boost::math::tools::evaluate_polynomial(coef, m - 0.875);
}
case 18:
case 19:
return boost::math::detail::ellint_k_imp(k, boost::math::policies::policy<>(), std::integral_constant<int, 2>());
default:
if (m == 1)
{
return boost::math::policies::raise_overflow_error<T>(function, nullptr, boost::math::policies::policy<>());
}
else
return boost::math::policies::raise_domain_error<T>(function,
"Got k = %1%, function requires |k| <= 1", k, boost::math::policies::policy<>());
}
}
template <typename Real>
void basic_ellint_rational_performance(benchmark::State& state)
{
std::vector<Real>const& test_set = get_test_data<Real>();
Real& r = get_result_data<Real>();
for(auto _ : state)
{
for(unsigned i = 0; i < test_set.size(); ++i)
benchmark::DoNotOptimize(r += basic_ellint_rational(test_set[i]));
}
}
template <typename Real>
void basic_ellint_rational_performance_no_error_check(benchmark::State& state)
{
std::vector<Real>const& test_set = get_test_data<Real>();
Real& r = get_result_data<Real>();
for(auto _ : state)
{
for(unsigned i = 0; i < test_set.size(); ++i)
benchmark::DoNotOptimize(r += basic_ellint_rational_no_error_checks(test_set[i]));
}
}
template <typename Real>
void ellint_1_performance(benchmark::State& state)
{
std::vector<Real>const& test_set = get_test_data<Real>();
Real& r = get_result_data<Real>();
for(auto _ : state)
{
for(unsigned i = 0; i < test_set.size(); ++i)
benchmark::DoNotOptimize(r += boost::math::ellint_1(test_set[i]));
}
}
template <typename Real>
void ellint_1_performance_no_overflow_check(benchmark::State& state)
{
std::vector<Real>const& test_set = get_test_data<Real>();
Real& r = get_result_data<Real>();
for(auto _ : state)
{
for(unsigned i = 0; i < test_set.size(); ++i)
benchmark::DoNotOptimize(r += boost::math::ellint_1(test_set[i], boost::math::policies::make_policy(boost::math::policies::overflow_error<boost::math::policies::ignore_error>())));
}
}
template <typename Real>
void ellint_2_performance(benchmark::State& state)
{
std::vector<Real>const& test_set = get_test_data<Real>();
Real& r = get_result_data<Real>();
for(auto _ : state)
{
for(unsigned i = 0; i < test_set.size(); ++i)
benchmark::DoNotOptimize(r += boost::math::ellint_2(test_set[i]));
}
}
template <typename Real>
void ellint_2_performance_no_overflow_check(benchmark::State& state)
{
std::vector<Real>const& test_set = get_test_data<Real>();
Real& r = get_result_data<Real>();
for(auto _ : state)
{
for(unsigned i = 0; i < test_set.size(); ++i)
benchmark::DoNotOptimize(r += boost::math::ellint_2(test_set[i], boost::math::policies::make_policy(boost::math::policies::overflow_error<boost::math::policies::ignore_error>())));
}
}
template <typename Real>
void log_performance(benchmark::State& state)
{
std::vector<Real>const& test_set = get_test_data<Real>();
Real& r = get_result_data<Real>();
for(auto _ : state)
{
for(unsigned i = 0; i < test_set.size(); ++i)
benchmark::DoNotOptimize(r += log(test_set[i]));
}
}
BENCHMARK_TEMPLATE(ellint_1_performance, float);
BENCHMARK_TEMPLATE(ellint_1_performance, double);
BENCHMARK_TEMPLATE(ellint_1_performance, long double);
BENCHMARK_TEMPLATE(ellint_1_performance_no_overflow_check, float);
BENCHMARK_TEMPLATE(ellint_1_performance_no_overflow_check, double);
BENCHMARK_TEMPLATE(ellint_1_performance_no_overflow_check, long double);
BENCHMARK_TEMPLATE(ellint_2_performance, float);
BENCHMARK_TEMPLATE(ellint_2_performance, double);
BENCHMARK_TEMPLATE(ellint_2_performance, long double);
BENCHMARK_TEMPLATE(ellint_2_performance_no_overflow_check, float);
BENCHMARK_TEMPLATE(ellint_2_performance_no_overflow_check, double);
BENCHMARK_TEMPLATE(ellint_2_performance_no_overflow_check, long double);
BENCHMARK_TEMPLATE(log_performance, float);
BENCHMARK_TEMPLATE(log_performance, double);
BENCHMARK_TEMPLATE(log_performance, long double);
BENCHMARK_TEMPLATE(basic_ellint_rational_performance, float);
BENCHMARK_TEMPLATE(basic_ellint_rational_performance, double);
BENCHMARK_TEMPLATE(basic_ellint_rational_performance, long double);
BENCHMARK_MAIN();