From 483eb2f56657e8e7f419ab1a4fab8dce9ade8609 Mon Sep 17 00:00:00 2001 From: Daniel Baumann Date: Sat, 27 Apr 2024 20:24:20 +0200 Subject: Adding upstream version 14.2.21. Signed-off-by: Daniel Baumann --- src/boost/libs/math/minimax/main.cpp | 650 +++++++++++++++++++++++++++++++++++ 1 file changed, 650 insertions(+) create mode 100644 src/boost/libs/math/minimax/main.cpp (limited to 'src/boost/libs/math/minimax/main.cpp') diff --git a/src/boost/libs/math/minimax/main.cpp b/src/boost/libs/math/minimax/main.cpp new file mode 100644 index 00000000..6ff01876 --- /dev/null +++ b/src/boost/libs/math/minimax/main.cpp @@ -0,0 +1,650 @@ +// (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_TEST_MODULE foobar +#define BOOST_UBLAS_TYPE_CHECK_EPSILON (type_traits::type_sqrt (boost::math::tools::epsilon ())) +#define BOOST_UBLAS_TYPE_CHECK_MIN (type_traits::type_sqrt ( boost::math::tools::min_value())) +#define BOOST_UBLAS_NDEBUG + +#include "multiprecision.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // for test_main +#include + + +extern mp_type f(const mp_type& x, int variant); +extern void show_extra( + const boost::math::tools::polynomial& n, + const boost::math::tools::polynomial& d, + const mp_type& x_offset, + const mp_type& y_offset, + int variant); + +using namespace boost::spirit::classic; + +mp_type 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(); +int working_precision = target_precision * 2; +bool started(false); +int variant(0); +int skew(0); +int brake(50); +mp_type x_offset(0), y_offset(0), x_scale(1); +bool auto_offset_y; + +boost::shared_ptr > p_remez; + +mp_type the_function(const mp_type& val) +{ + return f(x_scale * (val + x_offset), variant) + y_offset; +} + +void step_some(unsigned count) +{ + try{ + set_working_precision(working_precision); + if(!started) + { + // + // If we have an automatic y-offset calculate it now: + // + if(auto_offset_y) + { + mp_type 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; + set_output_precision(5); + std::cout << "Setting auto-y-offset to " << y_offset << std::endl; + } + // + // Truncate offsets to float precision: + // + x_offset = round_to_precision(x_offset, 20); + y_offset = round_to_precision(y_offset, 20); + // + // Construct new Remez state machine: + // + p_remez.reset(new boost::math::tools::remez_minimax( + &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(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); + mp_type r = p_remez->iterate(); + set_output_precision(3); + std::cout + << "Maximum Deviation Found: " << std::setprecision(3) << std::scientific << boost::math::tools::real_cast(p_remez->max_error()) << std::endl + << "Expected Error Term: " << std::setprecision(3) << std::scientific << boost::math::tools::real_cast(p_remez->error_term()) << std::endl + << "Maximum Relative Change in Control Points: " << std::setprecision(3) << std::scientific << boost::math::tools::real_cast(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*) +{ + set_working_precision(working_precision); + if(started) + { + boost::math::tools::polynomial n = p_remez->numerator(); + boost::math::tools::polynomial d = p_remez->denominator(); + std::vector cn = n.chebyshev(); + std::vector cd = d.chebyshev(); + int prec = 2 + (target_precision * 3010LL)/10000; + std::cout << std::scientific << std::setprecision(prec); + set_output_precision(prec); + boost::numeric::ublas::vector 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) +{ + set_working_precision(working_precision); + mp_type step = (b - a) / (points - 1); + mp_type x = a; + while(points > 1) + { + set_output_precision(10); + std::cout << std::setprecision(10) << std::setw(30) << std::left + << boost::lexical_cast(x) << the_function(x) << std::endl; + --points; + x += step; + } + std::cout << std::setprecision(10) << std::setw(30) << std::left + << boost::lexical_cast(b) << the_function(b) << std::endl; +} + +void graph(const char*, const char*) +{ + do_graph(3); +} + +template +mp_type convert_to_rr(const T& val) +{ + return val; +} +template +mp_type convert_to_rr(const boost::multiprecision::number& val) +{ + return boost::lexical_cast(val.str()); +} + +template +void do_test(T, const char* name) +{ + set_working_precision(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 n, d; + boost::math::tools::polynomial nr, dr; + nr = p_remez->numerator(); + dr = p_remez->denominator(); + n = nr; + d = dr; + + std::vector cn1, cd1; + cn1 = nr.chebyshev(); + cd1 = dr.chebyshev(); + std::vector cn, cd; + for(unsigned i = 0; i < cn1.size(); ++i) + { + cn.push_back(boost::math::tools::real_cast(cn1[i])); + } + for(unsigned i = 0; i < cd1.size(); ++i) + { + cd.push_back(boost::math::tools::real_cast(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 + zeros(p_remez->zero_points()), + cheb(p_remez->chebyshev_points()); + + mp_type 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) + { + mp_type true_result = the_function(zeros[i]); + T absissa = boost::math::tools::real_cast(zeros[i]); + mp_type test_result = convert_to_rr(n.evaluate(absissa) / d.evaluate(absissa)); + mp_type cheb_result = convert_to_rr(boost::math::tools::evaluate_chebyshev(cn, absissa) / boost::math::tools::evaluate_chebyshev(cd, absissa)); + mp_type 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(err) << boost::math::tools::real_cast(cheb_err) << std::endl; + } + // + // Do the tests at the Chebeshev control points: + // + for(unsigned i = 0; i < cheb.size(); ++i) + { + mp_type true_result = the_function(cheb[i]); + T absissa = boost::math::tools::real_cast(cheb[i]); + mp_type test_result = convert_to_rr(n.evaluate(absissa) / d.evaluate(absissa)); + mp_type cheb_result = convert_to_rr(boost::math::tools::evaluate_chebyshev(cn, absissa) / boost::math::tools::evaluate_chebyshev(cd, absissa)); + mp_type 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(err) << + boost::math::tools::real_cast(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(max_error) << "Cheb: " << boost::math::tools::real_cast(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_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"); + do_test(double(0), "double"); + do_test((long double)(0), "long double"); +} + +template +void do_test_n(T, const char* name, unsigned count) +{ + set_working_precision(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 n, d; + boost::math::tools::polynomial nr, dr; + nr = p_remez->numerator(); + dr = p_remez->denominator(); + n = nr; + d = dr; + + std::vector cn1, cd1; + cn1 = nr.chebyshev(); + cd1 = dr.chebyshev(); + std::vector cn, cd; + for(unsigned i = 0; i < cn1.size(); ++i) + { + cn.push_back(boost::math::tools::real_cast(cn1[i])); + } + for(unsigned i = 0; i < cd1.size(); ++i) + { + cd.push_back(boost::math::tools::real_cast(cd1[i])); + } + + mp_type max_error(0), max_cheb_error(0); + mp_type 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(mp_type x = a; x <= b; x += step) + { + mp_type true_result = the_function(x); + //std::cout << true_result << std::endl; + T absissa = boost::math::tools::real_cast(x); + mp_type test_result = convert_to_rr(n.evaluate(absissa) / d.evaluate(absissa)); + //std::cout << test_result << std::endl; + mp_type 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; + mp_type 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(absissa) + << (test_result < true_result ? "-" : "") << std::setw(20) << std::left + << boost::math::tools::real_cast(err) + << boost::math::tools::real_cast(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(max_error) + << " Cheb: " << boost::math::tools::real_cast(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(mp_type(), "mp_type", 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 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) + { + 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; + set_working_precision(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 n, d; + n = p_remez->numerator(); + d = p_remez->denominator(); + + mp_type max_error(0); + mp_type step = (b - a) / i; + + std::cout << "Evaluating Numerator...\n"; + mp_type 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; + } +} + +BOOST_AUTO_TEST_CASE( test_main ) +{ + std::string line; + real_parser const rr_p; + while(std::getline(std::cin, line)) + { + if(parse(line.c_str(), str_p("quit"), space_p).full) + return; + 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("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("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; + } + } +} -- cgit v1.2.3