summaryrefslogtreecommitdiff
path: root/libs/math/minimax
diff options
context:
space:
mode:
Diffstat (limited to 'libs/math/minimax')
-rw-r--r--libs/math/minimax/f.cpp49
-rw-r--r--libs/math/minimax/main.cpp82
2 files changed, 92 insertions, 39 deletions
diff --git a/libs/math/minimax/f.cpp b/libs/math/minimax/f.cpp
index 3435f0175..cc52b1fcc 100644
--- a/libs/math/minimax/f.cpp
+++ b/libs/math/minimax/f.cpp
@@ -194,42 +194,13 @@ boost::math::ntl::RR f(const boost::math::ntl::RR& x, int variant)
case 13:
// K(k):
{
- // x = k^2
- boost::math::ntl::RR k2 = x;
- if(k2 > boost::math::ntl::RR(1) - 1e-40)
- k2 = boost::math::ntl::RR(1) - 1e-40;
- /*if(k < 1e-40)
- k = 1e-40;*/
- boost::math::ntl::RR p2 = boost::math::constants::pi<boost::math::ntl::RR>() / 2;
- return (boost::math::ellint_1(sqrt(k2))) / (p2 - boost::math::log1p(-k2));
+ return boost::math::ellint_1(x);
}
case 14:
// K(k)
{
- static double P[] =
- {
- 1.38629436111989062502E0,
- 9.65735902811690126535E-2,
- 3.08851465246711995998E-2,
- 1.49380448916805252718E-2,
- 8.79078273952743772254E-3,
- 6.18901033637687613229E-3,
- 6.87489687449949877925E-3,
- 9.85821379021226008714E-3,
- 7.97404013220415179367E-3,
- 2.28025724005875567385E-3,
- 1.37982864606273237150E-4
- };
-
- // x = 1 - k^2
- boost::math::ntl::RR mp = x;
- if(mp < 1e-20)
- mp = 1e-20;
- if(mp == 1)
- mp -= 1e-20;
- boost::math::ntl::RR k = sqrt(1 - mp);
- return (boost::math::ellint_1(k) - mp/*boost::math::tools::evaluate_polynomial(P, mp)*/) / -log(mp);
- }
+ return boost::math::ellint_1(1-x) / log(x);
+ }
case 15:
// E(k)
{
@@ -323,6 +294,20 @@ boost::math::ntl::RR f(const boost::math::ntl::RR& x, int variant)
// cbrt over [0.5, 1]
return boost::math::cbrt(x);
}
+ case 30:
+ {
+ // trigamma over [x,y]
+ boost::math::ntl::RR y = x;
+ y = sqrt(y);
+ return boost::math::trigamma(x) * (x * x);
+ }
+ case 31:
+ {
+ // trigamma over [x, INF]
+ if(x == 0) return 1;
+ boost::math::ntl::RR y = (x == 0) ? (std::numeric_limits<double>::max)() / 2 : 1/x;
+ return boost::math::trigamma(y) * y;
+ }
}
return 0;
}
diff --git a/libs/math/minimax/main.cpp b/libs/math/minimax/main.cpp
index a12657ed8..bc430ec4f 100644
--- a/libs/math/minimax/main.cpp
+++ b/libs/math/minimax/main.cpp
@@ -3,6 +3,7 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#define BOOST_TEST_MODULE foobar
#define BOOST_UBLAS_TYPE_CHECK_EPSILON (type_traits<real_type>::type_sqrt (boost::math::tools::epsilon <real_type>()))
#define BOOST_UBLAS_TYPE_CHECK_MIN (type_traits<real_type>::type_sqrt ( boost::math::tools::min_value<real_type>()))
#define BOOST_UBLAS_NDEBUG
@@ -21,6 +22,7 @@ using boost::math::ntl::pow;
#include <iomanip>
#include <string>
#include <boost/test/included/unit_test.hpp> // for test_main
+#include <boost/multiprecision/cpp_bin_float.hpp>
extern boost::math::ntl::RR f(const boost::math::ntl::RR& x, int variant);
extern void show_extra(
@@ -211,6 +213,41 @@ void graph(const char*, const char*)
}
template <class T>
+boost::math::ntl::RR convert_to_rr(const T& val)
+{
+ return val;
+}
+template <class Backend, boost::multiprecision::expression_template_option ET>
+boost::math::ntl::RR convert_to_rr(const boost::multiprecision::number<Backend, ET>& val)
+{
+ return boost::lexical_cast<boost::math::ntl::RR>(val.str());
+}
+
+namespace boost{ namespace math{ namespace tools{
+
+template <>
+boost::multiprecision::cpp_bin_float_double_extended real_cast<boost::multiprecision::cpp_bin_float_double_extended, boost::math::ntl::RR>(boost::math::ntl::RR val)
+{
+ unsigned p = NTL::RR::OutputPrecision();
+ NTL::RR::SetOutputPrecision(20);
+ boost::multiprecision::cpp_bin_float_double_extended r = boost::lexical_cast<boost::multiprecision::cpp_bin_float_double_extended>(val);
+ NTL::RR::SetOutputPrecision(p);
+ return r;
+}
+template <>
+boost::multiprecision::cpp_bin_float_quad real_cast<boost::multiprecision::cpp_bin_float_quad, boost::math::ntl::RR>(boost::math::ntl::RR val)
+{
+ unsigned p = NTL::RR::OutputPrecision();
+ NTL::RR::SetOutputPrecision(35);
+ boost::multiprecision::cpp_bin_float_quad r = boost::lexical_cast<boost::multiprecision::cpp_bin_float_quad>(val);
+ NTL::RR::SetOutputPrecision(p);
+ return r;
+}
+
+}}}
+
+
+template <class T>
void do_test(T, const char* name)
{
boost::math::ntl::RR::SetPrecision(working_precision);
@@ -260,8 +297,8 @@ void do_test(T, const char* name)
{
boost::math::ntl::RR true_result = the_function(zeros[i]);
T absissa = boost::math::tools::real_cast<T>(zeros[i]);
- boost::math::ntl::RR test_result = n.evaluate(absissa) / d.evaluate(absissa);
- boost::math::ntl::RR cheb_result = boost::math::tools::evaluate_chebyshev(cn, absissa) / boost::math::tools::evaluate_chebyshev(cd, absissa);
+ boost::math::ntl::RR test_result = convert_to_rr(n.evaluate(absissa) / d.evaluate(absissa));
+ boost::math::ntl::RR cheb_result = convert_to_rr(boost::math::tools::evaluate_chebyshev(cn, absissa) / boost::math::tools::evaluate_chebyshev(cd, absissa));
boost::math::ntl::RR err, cheb_err;
if(rel_error)
{
@@ -287,8 +324,8 @@ void do_test(T, const char* name)
{
boost::math::ntl::RR true_result = the_function(cheb[i]);
T absissa = boost::math::tools::real_cast<T>(cheb[i]);
- boost::math::ntl::RR test_result = n.evaluate(absissa) / d.evaluate(absissa);
- boost::math::ntl::RR cheb_result = boost::math::tools::evaluate_chebyshev(cn, absissa) / boost::math::tools::evaluate_chebyshev(cd, absissa);
+ boost::math::ntl::RR test_result = convert_to_rr(n.evaluate(absissa) / d.evaluate(absissa));
+ boost::math::ntl::RR cheb_result = convert_to_rr(boost::math::tools::evaluate_chebyshev(cn, absissa) / boost::math::tools::evaluate_chebyshev(cd, absissa));
boost::math::ntl::RR err, cheb_err;
if(rel_error)
{
@@ -334,6 +371,16 @@ void test_long(const char*, const char*)
do_test((long double)(0), "long double");
}
+void test_float80(const char*, const char*)
+{
+ do_test((boost::multiprecision::cpp_bin_float_double_extended)(0), "float80");
+}
+
+void test_float128(const char*, const char*)
+{
+ do_test((boost::multiprecision::cpp_bin_float_quad)(0), "float128");
+}
+
void test_all(const char*, const char*)
{
do_test(float(0), "float");
@@ -383,9 +430,12 @@ void do_test_n(T, const char* name, unsigned count)
for(boost::math::ntl::RR x = a; x <= b; x += step)
{
boost::math::ntl::RR true_result = the_function(x);
+ //std::cout << true_result << std::endl;
T absissa = boost::math::tools::real_cast<T>(x);
- boost::math::ntl::RR test_result = n.evaluate(absissa) / d.evaluate(absissa);
- boost::math::ntl::RR cheb_result = boost::math::tools::evaluate_chebyshev(cn, absissa) / boost::math::tools::evaluate_chebyshev(cd, absissa);
+ boost::math::ntl::RR test_result = convert_to_rr(n.evaluate(absissa) / d.evaluate(absissa));
+ //std::cout << test_result << std::endl;
+ boost::math::ntl::RR cheb_result = convert_to_rr(boost::math::tools::evaluate_chebyshev(cn, absissa) / boost::math::tools::evaluate_chebyshev(cd, absissa));
+ //std::cout << cheb_result << std::endl;
boost::math::ntl::RR err, cheb_err;
if(rel_error)
{
@@ -441,6 +491,16 @@ void test_long_n(unsigned n)
do_test_n((long double)(0), "long double", n);
}
+void test_float80_n(unsigned n)
+{
+ do_test_n((boost::multiprecision::cpp_bin_float_double_extended)(0), "float80", n);
+}
+
+void test_float128_n(unsigned n)
+{
+ do_test_n((boost::multiprecision::cpp_bin_float_quad)(0), "float128", n);
+}
+
void rotate(const char*, const char*)
{
if(p_remez)
@@ -504,7 +564,7 @@ BOOST_AUTO_TEST_CASE( test_main )
while(std::getline(std::cin, line))
{
if(parse(line.c_str(), str_p("quit"), space_p).full)
- return 0;
+ return;
if(false == parse(line.c_str(),
(
@@ -570,6 +630,14 @@ BOOST_AUTO_TEST_CASE( test_main )
||
str_p("test") && str_p("long")[&test_long]
||
+ str_p("test") && str_p("float80") && uint_p[&test_float80_n]
+ ||
+ str_p("test") && str_p("float80")[&test_float80]
+ ||
+ str_p("test") && str_p("float128") && uint_p[&test_float128_n]
+ ||
+ str_p("test") && str_p("float128")[&test_float128]
+ ||
str_p("test") && str_p("all")[&test_all]
||
str_p("test") && uint_p[&test_n]
View Lorry Controller Status