| // (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) |
| |
| #include <pch.hpp> |
| |
| #include <boost/test/test_exec_monitor.hpp> |
| #include <boost/test/floating_point_comparison.hpp> |
| |
| #include <boost/math/tools/toms748_solve.hpp> |
| #include <boost/math/special_functions/gamma.hpp> |
| #include <boost/math/special_functions/beta.hpp> |
| |
| // |
| // Test functor implements the same test cases as used by |
| // "Algorithm 748: Enclosing Zeros of Continuous Functions" |
| // by G.E. Alefeld, F.A. Potra and Yixun Shi. |
| // |
| // Plus two more: one for inverting the incomplete gamma, |
| // and one for inverting the incomplete beta. |
| // |
| template <class T> |
| struct toms748tester |
| { |
| toms748tester(unsigned i) : k(i), ip(0), a(0), b(0){} |
| toms748tester(unsigned i, int ip_) : k(i), ip(ip_), a(0), b(0){} |
| toms748tester(unsigned i, T a_, T b_) : k(i), ip(0), a(a_), b(b_){} |
| |
| static unsigned total_calls() |
| { |
| return invocations; |
| } |
| static void reset() |
| { |
| invocations = 0; |
| } |
| |
| T operator()(T x) |
| { |
| using namespace std; |
| |
| ++invocations; |
| switch(k) |
| { |
| case 1: |
| return sin(x) - x / 2; |
| case 2: |
| { |
| T r = 0; |
| for(int i = 1; i <= 20; ++i) |
| { |
| T p = (2 * i - 5); |
| T q = x - i * i; |
| r += p * p / (q * q * q); |
| } |
| r *= -2; |
| return r; |
| } |
| case 3: |
| return a * x * exp(b * x); |
| case 4: |
| return pow(x, b) - a; |
| case 5: |
| return sin(x) - 0.5; |
| case 6: |
| return 2 * x * exp(-T(ip)) - 2 * exp(-ip * x) + 1; |
| case 7: |
| return (1 + (1 - ip) * (1 - ip)) * x - (1 - ip * x) * (1 - ip * x); |
| case 8: |
| return x * x - pow(1 - x, a); |
| case 9: |
| return (1 + (1 - ip) * (1 - ip) * (1 - ip) * (1 - ip)) * x - (1 - ip * x) * (1 - ip * x) * (1 - ip * x) * (1 - ip * x); |
| case 10: |
| return exp(-ip * x) * (x - 1) + pow(x, T(ip)); |
| case 11: |
| return (ip * x - 1) / ((ip - 1) * x); |
| case 12: |
| return pow(x, T(1)/ip) - pow(T(ip), T(1)/ip); |
| case 13: |
| return x == 0 ? 0 : x * exp(-1 / (x * x)); |
| case 14: |
| return x >= 0 ? (T(ip)/20) * (x / 1.5f + sin(x) - 1) : -T(ip)/20; |
| case 15: |
| { |
| T d = 2e-3f / (1+ip); |
| if(x > d) |
| return exp(1.0) - 1.859; |
| else if(x > 0) |
| return exp((ip+1)*x*1000 / 2) - 1.859; |
| return -0.859f; |
| } |
| case 16: |
| { |
| return boost::math::gamma_q(x, a) - b; |
| } |
| case 17: |
| return boost::math::ibeta(x, a, 0.5) - b; |
| } |
| return 0; |
| } |
| private: |
| int k; // index of problem. |
| int ip; // integer parameter |
| T a, b; // real parameter |
| |
| static unsigned invocations; |
| }; |
| |
| template <class T> |
| unsigned toms748tester<T>::invocations = 0; |
| |
| boost::uintmax_t total = 0; |
| boost::uintmax_t invocations = 0; |
| |
| template <class T> |
| void run_test(T a, T b, int id) |
| { |
| boost::uintmax_t c = 1000; |
| std::pair<double, double> r = toms748_solve(toms748tester<double>(id), |
| a, |
| b, |
| boost::math::tools::eps_tolerance<double>(std::numeric_limits<double>::digits), |
| c); |
| BOOST_CHECK_EQUAL(c, toms748tester<double>::total_calls()); |
| total += c; |
| invocations += toms748tester<double>::total_calls(); |
| std::cout << "Function " << id << "\nresult={" << r.first << ", " << r.second << "} total calls=" << toms748tester<double>::total_calls() << "\n\n"; |
| toms748tester<double>::reset(); |
| } |
| |
| template <class T> |
| void run_test(T a, T b, int id, int p) |
| { |
| boost::uintmax_t c = 1000; |
| std::pair<double, double> r = toms748_solve(toms748tester<double>(id, p), |
| a, |
| b, |
| boost::math::tools::eps_tolerance<double>(std::numeric_limits<double>::digits), |
| c); |
| BOOST_CHECK_EQUAL(c, toms748tester<double>::total_calls()); |
| total += c; |
| invocations += toms748tester<double>::total_calls(); |
| std::cout << "Function " << id << "\nresult={" << r.first << ", " << r.second << "} total calls=" << toms748tester<double>::total_calls() << "\n\n"; |
| toms748tester<double>::reset(); |
| } |
| |
| template <class T> |
| void run_test(T a, T b, int id, T p1, T p2) |
| { |
| boost::uintmax_t c = 1000; |
| std::pair<double, double> r = toms748_solve(toms748tester<double>(id, p1, p2), |
| a, |
| b, |
| boost::math::tools::eps_tolerance<double>(std::numeric_limits<double>::digits), |
| c); |
| BOOST_CHECK_EQUAL(c, toms748tester<double>::total_calls()); |
| total += c; |
| invocations += toms748tester<double>::total_calls(); |
| std::cout << "Function " << id << "\n Result={" << r.first << ", " << r.second << "} total calls=" << toms748tester<double>::total_calls() << "\n\n"; |
| toms748tester<double>::reset(); |
| } |
| |
| int test_main(int, char* []) |
| { |
| std::cout << std::setprecision(18); |
| run_test(3.14/2, 3.14, 1); |
| |
| for(int i = 1; i <= 10; i += 1) |
| { |
| run_test(i*i + 1e-9, (i+1)*(i+1) - 1e-9, 2); |
| } |
| |
| run_test(-9.0, 31.0, 3, -40.0, -1.0); |
| run_test(-9.0, 31.0, 3, -100.0, -2.0); |
| run_test(-9.0, 31.0, 3, -200.0, -3.0); |
| |
| for(int n = 4; n <= 12; n += 2) |
| { |
| run_test(0.0, 5.0, 4, 0.2, double(n)); |
| } |
| for(int n = 4; n <= 12; n += 2) |
| { |
| run_test(0.0, 5.0, 4, 1.0, double(n)); |
| } |
| for(int n = 8; n <= 14; n += 2) |
| { |
| run_test(-0.95, 4.05, 4, 1.0, double(n)); |
| } |
| run_test(0.0, 1.5, 5); |
| for(int n = 1; n <= 5; ++n) |
| { |
| run_test(0.0, 1.0, 6, n); |
| } |
| for(int n = 20; n <= 100; n += 20) |
| { |
| run_test(0.0, 1.0, 6, n); |
| } |
| run_test(0.0, 1.0, 7, 5); |
| run_test(0.0, 1.0, 7, 10); |
| run_test(0.0, 1.0, 7, 20); |
| run_test(0.0, 1.0, 8, 2); |
| run_test(0.0, 1.0, 8, 5); |
| run_test(0.0, 1.0, 8, 10); |
| run_test(0.0, 1.0, 8, 15); |
| run_test(0.0, 1.0, 8, 20); |
| run_test(0.0, 1.0, 9, 1); |
| run_test(0.0, 1.0, 9, 2); |
| run_test(0.0, 1.0, 9, 4); |
| run_test(0.0, 1.0, 9, 5); |
| run_test(0.0, 1.0, 9, 8); |
| run_test(0.0, 1.0, 9, 15); |
| run_test(0.0, 1.0, 9, 20); |
| run_test(0.0, 1.0, 10, 1); |
| run_test(0.0, 1.0, 10, 5); |
| run_test(0.0, 1.0, 10, 10); |
| run_test(0.0, 1.0, 10, 15); |
| run_test(0.0, 1.0, 10, 20); |
| |
| run_test(0.01, 1.0, 11, 2); |
| run_test(0.01, 1.0, 11, 5); |
| run_test(0.01, 1.0, 11, 15); |
| run_test(0.01, 1.0, 11, 20); |
| |
| for(int n = 2; n <= 6; ++n) |
| run_test(1.0, 100.0, 12, n); |
| for(int n = 7; n <= 33; n+=2) |
| run_test(1.0, 100.0, 12, n); |
| |
| run_test(-1.0, 4.0, 13); |
| |
| for(int n = 1; n <= 40; ++n) |
| run_test(-1e4, 3.14/2, 14, n); |
| |
| for(int n = 20; n <= 40; ++n) |
| run_test(-1e4, 1e-4, 15, n); |
| |
| for(int n = 100; n <= 1000; n+=100) |
| run_test(-1e4, 1e-4, 15, n); |
| |
| std::cout << "Total iterations consumed = " << total << std::endl; |
| std::cout << "Total function invocations consumed = " << invocations << std::endl << std::endl; |
| |
| BOOST_CHECK(invocations < 3150); |
| |
| std::cout << std::setprecision(18); |
| |
| for(int n = 5; n <= 100; n += 10) |
| run_test(sqrt(double(n)), double(n+1), 16, (double)n, 0.4); |
| |
| for(int n = 5; n <= 100; n += 10) |
| run_test(double(n / 2), double(2*n), 17, (double)n, 0.4); |
| |
| |
| for(int n = 4; n <= 12; n += 2) |
| { |
| boost::uintmax_t c = 1000; |
| std::pair<double, double> r = bracket_and_solve_root(toms748tester<double>(4, 0.2, double(n)), |
| 2.0, |
| 2.0, |
| true, |
| boost::math::tools::eps_tolerance<double>(std::numeric_limits<double>::digits), |
| c); |
| std::cout << std::setprecision(18); |
| std::cout << "Function " << 4 << "\n Result={" << r.first << ", " << r.second << "} total calls=" << toms748tester<double>::total_calls() << "\n\n"; |
| toms748tester<double>::reset(); |
| BOOST_CHECK(c < 20); |
| } |
| |
| return 0; |
| } |
| |