| // (C) Copyright John Maddock 2006. |
| // 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) |
| |
| #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 |
| |
| #include <boost/math/bindings/rr.hpp> |
| namespace std{ |
| using boost::math::ntl::pow; |
| } // workaround for spirit parser. |
| #include <boost/math/tools/remez.hpp> |
| #include <boost/math/tools/test.hpp> |
| #include <boost/math/special_functions/binomial.hpp> |
| #include <boost/spirit/core.hpp> |
| #include <boost/spirit/actor.hpp> |
| #include <boost/lexical_cast.hpp> |
| #include <iostream> |
| #include <iomanip> |
| #include <string> |
| #include <boost/test/included/test_exec_monitor.hpp> // for test_main |
| |
| extern boost::math::ntl::RR f(const boost::math::ntl::RR& x, int variant); |
| extern void show_extra( |
| const boost::math::tools::polynomial<boost::math::ntl::RR>& n, |
| const boost::math::tools::polynomial<boost::math::ntl::RR>& d, |
| const boost::math::ntl::RR& x_offset, |
| const boost::math::ntl::RR& y_offset, |
| int variant); |
| |
| using namespace boost::spirit; |
| |
| boost::math::ntl::RR a(0), b(1); // range to optimise over |
| bool rel_error(true); |
| bool pin(false); |
| int orderN(3); |
| int orderD(1); |
| int target_precision = boost::math::tools::digits<long double>(); |
| int working_precision = target_precision * 2; |
| bool started(false); |
| int variant(0); |
| int skew(0); |
| int brake(50); |
| boost::math::ntl::RR x_offset(0), y_offset(0), x_scale(1); |
| bool auto_offset_y; |
| |
| boost::shared_ptr<boost::math::tools::remez_minimax<boost::math::ntl::RR> > p_remez; |
| |
| boost::math::ntl::RR the_function(const boost::math::ntl::RR& val) |
| { |
| return f(x_scale * (val + x_offset), variant) + y_offset; |
| } |
| |
| void step_some(unsigned count) |
| { |
| try{ |
| NTL::RR::SetPrecision(working_precision); |
| if(!started) |
| { |
| // |
| // If we have an automatic y-offset calculate it now: |
| // |
| if(auto_offset_y) |
| { |
| boost::math::ntl::RR fa, fb, fm; |
| fa = f(x_scale * (a + x_offset), variant); |
| fb = f(x_scale * (b + x_offset), variant); |
| fm = f(x_scale * ((a+b)/2 + x_offset), variant); |
| y_offset = -(fa + fb + fm) / 3; |
| NTL::RR::SetOutputPrecision(5); |
| std::cout << "Setting auto-y-offset to " << y_offset << std::endl; |
| } |
| // |
| // Truncate offsets to float precision: |
| // |
| x_offset = NTL::RoundToPrecision(x_offset.value(), 20); |
| y_offset = NTL::RoundToPrecision(y_offset.value(), 20); |
| // |
| // Construct new Remez state machine: |
| // |
| p_remez.reset(new boost::math::tools::remez_minimax<boost::math::ntl::RR>( |
| &the_function, |
| orderN, orderD, |
| a, b, |
| pin, |
| rel_error, |
| skew, |
| working_precision)); |
| std::cout << "Max error in interpolated form: " << std::setprecision(3) << std::scientific << boost::math::tools::real_cast<double>(p_remez->max_error()) << std::endl; |
| // |
| // Signal that we've started: |
| // |
| started = true; |
| } |
| unsigned i; |
| for(i = 0; i < count; ++i) |
| { |
| std::cout << "Stepping..." << std::endl; |
| p_remez->set_brake(brake); |
| boost::math::ntl::RR r = p_remez->iterate(); |
| NTL::RR::SetOutputPrecision(3); |
| std::cout |
| << "Maximum Deviation Found: " << std::setprecision(3) << std::scientific << boost::math::tools::real_cast<double>(p_remez->max_error()) << std::endl |
| << "Expected Error Term: " << std::setprecision(3) << std::scientific << boost::math::tools::real_cast<double>(p_remez->error_term()) << std::endl |
| << "Maximum Relative Change in Control Points: " << std::setprecision(3) << std::scientific << boost::math::tools::real_cast<double>(r) << std::endl; |
| } |
| } |
| catch(const std::exception& e) |
| { |
| std::cout << "Step failed with exception: " << e.what() << std::endl; |
| } |
| } |
| |
| void step(const char*, const char*) |
| { |
| step_some(1); |
| } |
| |
| void show(const char*, const char*) |
| { |
| NTL::RR::SetPrecision(working_precision); |
| if(started) |
| { |
| boost::math::tools::polynomial<boost::math::ntl::RR> n = p_remez->numerator(); |
| boost::math::tools::polynomial<boost::math::ntl::RR> d = p_remez->denominator(); |
| std::vector<boost::math::ntl::RR> cn = n.chebyshev(); |
| std::vector<boost::math::ntl::RR> cd = d.chebyshev(); |
| int prec = 2 + (target_precision * 3010LL)/10000; |
| std::cout << std::scientific << std::setprecision(prec); |
| NTL::RR::SetOutputPrecision(prec); |
| boost::numeric::ublas::vector<boost::math::ntl::RR> v = p_remez->zero_points(); |
| |
| std::cout << " Zeros = {\n"; |
| unsigned i; |
| for(i = 0; i < v.size(); ++i) |
| { |
| std::cout << " " << v[i] << std::endl; |
| } |
| std::cout << " }\n"; |
| |
| v = p_remez->chebyshev_points(); |
| std::cout << " Chebeshev Control Points = {\n"; |
| for(i = 0; i < v.size(); ++i) |
| { |
| std::cout << " " << v[i] << std::endl; |
| } |
| std::cout << " }\n"; |
| |
| std::cout << "X offset: " << x_offset << std::endl; |
| std::cout << "X scale: " << x_scale << std::endl; |
| std::cout << "Y offset: " << y_offset << std::endl; |
| |
| std::cout << "P = {"; |
| for(i = 0; i < n.size(); ++i) |
| { |
| std::cout << " " << n[i] << "L," << std::endl; |
| } |
| std::cout << " }\n"; |
| |
| std::cout << "Q = {"; |
| for(i = 0; i < d.size(); ++i) |
| { |
| std::cout << " " << d[i] << "L," << std::endl; |
| } |
| std::cout << " }\n"; |
| |
| std::cout << "CP = {"; |
| for(i = 0; i < cn.size(); ++i) |
| { |
| std::cout << " " << cn[i] << "L," << std::endl; |
| } |
| std::cout << " }\n"; |
| |
| std::cout << "CQ = {"; |
| for(i = 0; i < cd.size(); ++i) |
| { |
| std::cout << " " << cd[i] << "L," << std::endl; |
| } |
| std::cout << " }\n"; |
| |
| show_extra(n, d, x_offset, y_offset, variant); |
| } |
| else |
| { |
| std::cerr << "Nothing to display" << std::endl; |
| } |
| } |
| |
| void do_graph(unsigned points) |
| { |
| NTL::RR::SetPrecision(working_precision); |
| boost::math::ntl::RR step = (b - a) / (points - 1); |
| boost::math::ntl::RR x = a; |
| while(points > 1) |
| { |
| NTL::RR::SetOutputPrecision(10); |
| std::cout << std::setprecision(10) << std::setw(30) << std::left |
| << boost::lexical_cast<std::string>(x) << the_function(x) << std::endl; |
| --points; |
| x += step; |
| } |
| std::cout << std::setprecision(10) << std::setw(30) << std::left |
| << boost::lexical_cast<std::string>(b) << the_function(b) << std::endl; |
| } |
| |
| void graph(const char*, const char*) |
| { |
| do_graph(3); |
| } |
| |
| template <class T> |
| void do_test(T, const char* name) |
| { |
| boost::math::ntl::RR::SetPrecision(working_precision); |
| if(started) |
| { |
| // |
| // We want to test the approximation at fixed precision: |
| // either float, double or long double. Begin by getting the |
| // polynomials: |
| // |
| boost::math::tools::polynomial<T> n, d; |
| boost::math::tools::polynomial<boost::math::ntl::RR> nr, dr; |
| nr = p_remez->numerator(); |
| dr = p_remez->denominator(); |
| n = nr; |
| d = dr; |
| |
| std::vector<boost::math::ntl::RR> cn1, cd1; |
| cn1 = nr.chebyshev(); |
| cd1 = dr.chebyshev(); |
| std::vector<T> cn, cd; |
| for(unsigned i = 0; i < cn1.size(); ++i) |
| { |
| cn.push_back(boost::math::tools::real_cast<T>(cn1[i])); |
| } |
| for(unsigned i = 0; i < cd1.size(); ++i) |
| { |
| cd.push_back(boost::math::tools::real_cast<T>(cd1[i])); |
| } |
| // |
| // We'll test at the Chebeshev control points which is where |
| // (in theory) the largest deviation should occur. For good |
| // measure we'll test at the zeros as well: |
| // |
| boost::numeric::ublas::vector<boost::math::ntl::RR> |
| zeros(p_remez->zero_points()), |
| cheb(p_remez->chebyshev_points()); |
| |
| boost::math::ntl::RR max_error(0), cheb_max_error(0); |
| |
| // |
| // Do the tests at the zeros: |
| // |
| std::cout << "Starting tests at " << name << " precision...\n"; |
| std::cout << "Absissa Error (Poly) Error (Cheb)\n"; |
| for(unsigned i = 0; i < zeros.size(); ++i) |
| { |
| 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 err, cheb_err; |
| if(rel_error) |
| { |
| err = boost::math::tools::relative_error(test_result, true_result); |
| cheb_err = boost::math::tools::relative_error(cheb_result, true_result); |
| } |
| else |
| { |
| err = fabs(test_result - true_result); |
| cheb_err = fabs(cheb_result - true_result); |
| } |
| if(err > max_error) |
| max_error = err; |
| if(cheb_err > cheb_max_error) |
| cheb_max_error = cheb_err; |
| std::cout << std::setprecision(6) << std::setw(15) << std::left << absissa |
| << std::setw(15) << std::left << boost::math::tools::real_cast<T>(err) << boost::math::tools::real_cast<T>(cheb_err) << std::endl; |
| } |
| // |
| // Do the tests at the Chebeshev control points: |
| // |
| for(unsigned i = 0; i < cheb.size(); ++i) |
| { |
| 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 err, cheb_err; |
| if(rel_error) |
| { |
| err = boost::math::tools::relative_error(test_result, true_result); |
| cheb_err = boost::math::tools::relative_error(cheb_result, true_result); |
| } |
| else |
| { |
| err = fabs(test_result - true_result); |
| cheb_err = fabs(cheb_result - true_result); |
| } |
| if(err > max_error) |
| max_error = err; |
| std::cout << std::setprecision(6) << std::setw(15) << std::left << absissa |
| << std::setw(15) << std::left << boost::math::tools::real_cast<T>(err) << |
| boost::math::tools::real_cast<T>(cheb_err) << std::endl; |
| } |
| std::string msg = "Max Error found at "; |
| msg += name; |
| msg += " precision = "; |
| msg.append(62 - 17 - msg.size(), ' '); |
| std::cout << msg << std::setprecision(6) << "Poly: " << std::setw(20) << std::left |
| << boost::math::tools::real_cast<T>(max_error) << "Cheb: " << boost::math::tools::real_cast<T>(cheb_max_error) << std::endl; |
| } |
| else |
| { |
| std::cout << "Nothing to test: try converging an approximation first!!!" << std::endl; |
| } |
| } |
| |
| void test_float(const char*, const char*) |
| { |
| do_test(float(0), "float"); |
| } |
| |
| void test_double(const char*, const char*) |
| { |
| do_test(double(0), "double"); |
| } |
| |
| void test_long(const char*, const char*) |
| { |
| do_test((long double)(0), "long double"); |
| } |
| |
| void test_all(const char*, const char*) |
| { |
| do_test(float(0), "float"); |
| do_test(double(0), "double"); |
| do_test((long double)(0), "long double"); |
| } |
| |
| template <class T> |
| void do_test_n(T, const char* name, unsigned count) |
| { |
| boost::math::ntl::RR::SetPrecision(working_precision); |
| if(started) |
| { |
| // |
| // We want to test the approximation at fixed precision: |
| // either float, double or long double. Begin by getting the |
| // polynomials: |
| // |
| boost::math::tools::polynomial<T> n, d; |
| boost::math::tools::polynomial<boost::math::ntl::RR> nr, dr; |
| nr = p_remez->numerator(); |
| dr = p_remez->denominator(); |
| n = nr; |
| d = dr; |
| |
| std::vector<boost::math::ntl::RR> cn1, cd1; |
| cn1 = nr.chebyshev(); |
| cd1 = dr.chebyshev(); |
| std::vector<T> cn, cd; |
| for(unsigned i = 0; i < cn1.size(); ++i) |
| { |
| cn.push_back(boost::math::tools::real_cast<T>(cn1[i])); |
| } |
| for(unsigned i = 0; i < cd1.size(); ++i) |
| { |
| cd.push_back(boost::math::tools::real_cast<T>(cd1[i])); |
| } |
| |
| boost::math::ntl::RR max_error(0), max_cheb_error(0); |
| boost::math::ntl::RR step = (b - a) / count; |
| |
| // |
| // Do the tests at the zeros: |
| // |
| std::cout << "Starting tests at " << name << " precision...\n"; |
| std::cout << "Absissa Error (poly) Error (Cheb)\n"; |
| for(boost::math::ntl::RR x = a; x <= b; x += step) |
| { |
| boost::math::ntl::RR true_result = the_function(x); |
| 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 err, cheb_err; |
| if(rel_error) |
| { |
| err = boost::math::tools::relative_error(test_result, true_result); |
| cheb_err = boost::math::tools::relative_error(cheb_result, true_result); |
| } |
| else |
| { |
| err = fabs(test_result - true_result); |
| cheb_err = fabs(cheb_result - true_result); |
| } |
| if(err > max_error) |
| max_error = err; |
| if(cheb_err > max_cheb_error) |
| max_cheb_error = cheb_err; |
| std::cout << std::setprecision(6) << std::setw(15) << std::left << boost::math::tools::real_cast<double>(absissa) |
| << (test_result < true_result ? "-" : "") << std::setw(20) << std::left |
| << boost::math::tools::real_cast<double>(err) |
| << boost::math::tools::real_cast<double>(cheb_err) << std::endl; |
| } |
| std::string msg = "Max Error found at "; |
| msg += name; |
| msg += " precision = "; |
| //msg.append(62 - 17 - msg.size(), ' '); |
| std::cout << msg << "Poly: " << std::setprecision(6) |
| //<< std::setw(15) << std::left |
| << boost::math::tools::real_cast<T>(max_error) |
| << " Cheb: " << boost::math::tools::real_cast<T>(max_cheb_error) << std::endl; |
| } |
| else |
| { |
| std::cout << "Nothing to test: try converging an approximation first!!!" << std::endl; |
| } |
| } |
| |
| void test_n(unsigned n) |
| { |
| do_test_n(boost::math::ntl::RR(), "boost::math::ntl::RR", n); |
| } |
| |
| void test_float_n(unsigned n) |
| { |
| do_test_n(float(0), "float", n); |
| } |
| |
| void test_double_n(unsigned n) |
| { |
| do_test_n(double(0), "double", n); |
| } |
| |
| void test_long_n(unsigned n) |
| { |
| do_test_n((long double)(0), "long double", n); |
| } |
| |
| void rotate(const char*, const char*) |
| { |
| if(p_remez) |
| { |
| p_remez->rotate(); |
| } |
| else |
| { |
| std::cerr << "Nothing to rotate" << std::endl; |
| } |
| } |
| |
| void rescale(const char*, const char*) |
| { |
| if(p_remez) |
| { |
| p_remez->rescale(a, b); |
| } |
| else |
| { |
| std::cerr << "Nothing to rescale" << std::endl; |
| } |
| } |
| |
| void graph_poly(const char*, const char*) |
| { |
| int i = 50; |
| boost::math::ntl::RR::SetPrecision(working_precision); |
| if(started) |
| { |
| // |
| // We want to test the approximation at fixed precision: |
| // either float, double or long double. Begin by getting the |
| // polynomials: |
| // |
| boost::math::tools::polynomial<boost::math::ntl::RR> n, d; |
| n = p_remez->numerator(); |
| d = p_remez->denominator(); |
| |
| boost::math::ntl::RR max_error(0); |
| boost::math::ntl::RR step = (b - a) / i; |
| |
| std::cout << "Evaluating Numerator...\n"; |
| boost::math::ntl::RR val; |
| for(val = a; val <= b; val += step) |
| std::cout << n.evaluate(val) << std::endl; |
| std::cout << "Evaluating Denominator...\n"; |
| for(val = a; val <= b; val += step) |
| std::cout << d.evaluate(val) << std::endl; |
| } |
| else |
| { |
| std::cout << "Nothing to test: try converging an approximation first!!!" << std::endl; |
| } |
| } |
| |
| int test_main(int, char* []) |
| { |
| std::string line; |
| real_parser<long double/*boost::math::ntl::RR*/ > const rr_p; |
| while(std::getline(std::cin, line)) |
| { |
| if(parse(line.c_str(), str_p("quit"), space_p).full) |
| return 0; |
| if(false == parse(line.c_str(), |
| ( |
| |
| str_p("range")[assign_a(started, false)] && real_p[assign_a(a)] && real_p[assign_a(b)] |
| || |
| str_p("relative")[assign_a(started, false)][assign_a(rel_error, true)] |
| || |
| str_p("absolute")[assign_a(started, false)][assign_a(rel_error, false)] |
| || |
| str_p("pin")[assign_a(started, false)] && str_p("true")[assign_a(pin, true)] |
| || |
| str_p("pin")[assign_a(started, false)] && str_p("false")[assign_a(pin, false)] |
| || |
| str_p("pin")[assign_a(started, false)] && str_p("1")[assign_a(pin, true)] |
| || |
| str_p("pin")[assign_a(started, false)] && str_p("0")[assign_a(pin, false)] |
| || |
| str_p("pin")[assign_a(started, false)][assign_a(pin, true)] |
| || |
| str_p("order")[assign_a(started, false)] && uint_p[assign_a(orderN)] && uint_p[assign_a(orderD)] |
| || |
| str_p("order")[assign_a(started, false)] && uint_p[assign_a(orderN)] |
| || |
| str_p("target-precision") && uint_p[assign_a(target_precision)] |
| || |
| str_p("working-precision")[assign_a(started, false)] && uint_p[assign_a(working_precision)] |
| || |
| str_p("variant")[assign_a(started, false)] && int_p[assign_a(variant)] |
| || |
| str_p("skew")[assign_a(started, false)] && int_p[assign_a(skew)] |
| || |
| str_p("brake") && int_p[assign_a(brake)] |
| || |
| str_p("step") && int_p[&step_some] |
| || |
| str_p("step")[&step] |
| || |
| str_p("poly")[&graph_poly] |
| || |
| str_p("info")[&show] |
| || |
| str_p("graph") && uint_p[&do_graph] |
| || |
| str_p("graph")[&graph] |
| || |
| str_p("x-offset") && real_p[assign_a(x_offset)] |
| || |
| str_p("x-scale") && real_p[assign_a(x_scale)] |
| || |
| str_p("y-offset") && str_p("auto")[assign_a(auto_offset_y, true)] |
| || |
| str_p("y-offset") && real_p[assign_a(y_offset)][assign_a(auto_offset_y, false)] |
| || |
| str_p("test") && str_p("float") && uint_p[&test_float_n] |
| || |
| str_p("test") && str_p("float")[&test_float] |
| || |
| str_p("test") && str_p("double") && uint_p[&test_double_n] |
| || |
| str_p("test") && str_p("double")[&test_double] |
| || |
| str_p("test") && str_p("long") && uint_p[&test_long_n] |
| || |
| str_p("test") && str_p("long")[&test_long] |
| || |
| str_p("test") && str_p("all")[&test_all] |
| || |
| str_p("test") && uint_p[&test_n] |
| || |
| str_p("rotate")[&rotate] |
| || |
| str_p("rescale") && real_p[assign_a(a)] && real_p[assign_a(b)] && epsilon_p[&rescale] |
| |
| ), space_p).full) |
| { |
| std::cout << "Unable to parse directive: \"" << line << "\"" << std::endl; |
| } |
| else |
| { |
| std::cout << "Variant = " << variant << std::endl; |
| std::cout << "range = [" << a << "," << b << "]" << std::endl; |
| std::cout << "Relative Error = " << rel_error << std::endl; |
| std::cout << "Pin to Origin = " << pin << std::endl; |
| std::cout << "Order (Num/Denom) = " << orderN << "/" << orderD << std::endl; |
| std::cout << "Target Precision = " << target_precision << std::endl; |
| std::cout << "Working Precision = " << working_precision << std::endl; |
| std::cout << "Skew = " << skew << std::endl; |
| std::cout << "Brake = " << brake << std::endl; |
| std::cout << "X Offset = " << x_offset << std::endl; |
| std::cout << "X scale = " << x_scale << std::endl; |
| std::cout << "Y Offset = "; |
| if(auto_offset_y) |
| std::cout << "Auto ("; |
| std::cout << y_offset; |
| if(auto_offset_y) |
| std::cout << ")"; |
| std::cout << std::endl; |
| } |
| } |
| return 0; |
| } |
