libs/math/test/test_toms748_solve.cpp

   1 //  (C) Copyright John Maddock 2006.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #include <pch.hpp>
   7
   8 #define BOOST_TEST_MAIN
   9 #include <boost/test/unit_test.hpp>
  10 #include <boost/test/tools/floating_point_comparison.hpp>
  11
  12 #include <boost/math/tools/toms748_solve.hpp>
  13 #include <boost/math/special_functions/gamma.hpp>
  14 #include <boost/math/special_functions/beta.hpp>
  15 #include <iostream>
  16 #include <iomanip>
  17
  18 //
  19 // Test functor implements the same test cases as used by
  20 // "Algorithm 748: Enclosing Zeros of Continuous Functions"
  21 // by G.E. Alefeld, F.A. Potra and Yixun Shi.
  22 //
  23 // Plus two more: one for inverting the incomplete gamma,
  24 // and one for inverting the incomplete beta.
  25 //
  26 template <class T>
  27 struct toms748tester
  28 {
  29    toms748tester(unsigned i) : k(i), ip(0), a(0), b(0){}
  30    toms748tester(unsigned i, int ip_) : k(i), ip(ip_), a(0), b(0){}
  31    toms748tester(unsigned i, T a_, T b_) : k(i), ip(0), a(a_), b(b_){}
  32
  33    static unsigned total_calls()
  34    {
  35       return invocations;
  36    }
  37    static void reset()
  38    {
  39       invocations = 0;
  40    }
  41
  42    T operator()(T x)
  43    {
  44       using namespace std;
  45
  46       ++invocations;
  47       switch(k)
  48       {
  49       case 1:
  50          return sin(x) - x / 2;
  51       case 2:
  52          {
  53             T r = 0;
  54             for(int i = 1; i <= 20; ++i)
  55             {
  56                T p = (2 * i - 5);
  57                T q = x - i * i;
  58                r += p * p / (q * q * q);
  59             }
  60             r *= -2;
  61             return r;
  62          }
  63       case 3:
  64          return a * x * exp(b * x);
  65       case 4:
  66          return pow(x, b) - a;
  67       case 5:
  68          return sin(x) - 0.5;
  69       case 6:
  70          return 2 * x * exp(-T(ip)) - 2 * exp(-ip * x) + 1;
  71       case 7:
  72          return (1 + (1 - ip) * (1 - ip)) * x - (1 - ip * x) * (1 - ip * x);
  73       case 8:
  74          return x * x - pow(1 - x, a);
  75       case 9:
  76          return (1 + (1 - ip) * (1 - ip) * (1 - ip) * (1 - ip)) * x - (1 - ip * x) * (1 - ip * x) * (1 - ip * x) * (1 - ip * x);
  77       case 10:
  78          return exp(-ip * x) * (x - 1) + pow(x, T(ip));
  79       case 11:
  80          return (ip * x - 1) / ((ip - 1) * x);
  81       case 12:
  82          return pow(x, T(1)/ip) - pow(T(ip), T(1)/ip);
  83       case 13:
  84          return x == 0 ? 0 : x * exp(-1 / (x * x));
  85       case 14:
  86          return x >= 0 ? (T(ip)/20) * (x / 1.5f + sin(x) - 1) : -T(ip)/20;
  87       case 15:
  88          {
  89             T d = 2e-3f / (1+ip);
  90             if(x > d)
  91                return exp(1.0) - 1.859;
  92             else if(x > 0)
  93                return exp((ip+1)*x*1000 / 2) - 1.859;
  94             return -0.859f;
  95          }
  96       case 16:
  97          {
  98             return boost::math::gamma_q(x, a) - b;
  99          }
 100       case 17:
 101          return boost::math::ibeta(x, a, 0.5) - b;
 102       }
 103       return 0;
 104    }
 105 private:
 106    int k; // index of problem.
 107    int ip; // integer parameter
 108    T a, b; // real parameter
 109
 110    static unsigned invocations;
 111 };
 112
 113 template <class T>
 114 unsigned toms748tester<T>::invocations = 0;
 115
 116 boost::uintmax_t total = 0;
 117 boost::uintmax_t invocations = 0;
 118
 119 template <class T>
 120 void run_test(T a, T b, int id)
 121 {
 122    boost::uintmax_t c = 1000;
 123    std::pair<double, double> r = toms748_solve(toms748tester<double>(id),
 124       a,
 125       b,
 126       boost::math::tools::eps_tolerance<double>(std::numeric_limits<double>::digits),
 127       c);
 128    BOOST_CHECK_EQUAL(c, toms748tester<double>::total_calls());
 129    total += c;
 130    invocations += toms748tester<double>::total_calls();
 131    std::cout << "Function " << id << "\nresult={" << r.first << ", " << r.second << "} total calls=" << toms748tester<double>::total_calls() << "\n\n";
 132    toms748tester<double>::reset();
 133 }
 134
 135 template <class T>
 136 void run_test(T a, T b, int id, int p)
 137 {
 138    boost::uintmax_t c = 1000;
 139    std::pair<double, double> r = toms748_solve(toms748tester<double>(id, p),
 140       a,
 141       b,
 142       boost::math::tools::eps_tolerance<double>(std::numeric_limits<double>::digits),
 143       c);
 144    BOOST_CHECK_EQUAL(c, toms748tester<double>::total_calls());
 145    total += c;
 146    invocations += toms748tester<double>::total_calls();
 147    std::cout << "Function " << id << "\nresult={" << r.first << ", " << r.second << "} total calls=" << toms748tester<double>::total_calls() << "\n\n";
 148    toms748tester<double>::reset();
 149 }
 150
 151 template <class T>
 152 void run_test(T a, T b, int id, T p1, T p2)
 153 {
 154    boost::uintmax_t c = 1000;
 155    std::pair<double, double> r = toms748_solve(toms748tester<double>(id, p1, p2),
 156       a,
 157       b,
 158       boost::math::tools::eps_tolerance<double>(std::numeric_limits<double>::digits),
 159       c);
 160    BOOST_CHECK_EQUAL(c, toms748tester<double>::total_calls());
 161    total += c;
 162    invocations += toms748tester<double>::total_calls();
 163    std::cout << "Function " << id << "\n   Result={" << r.first << ", " << r.second << "} total calls=" << toms748tester<double>::total_calls() << "\n\n";
 164    toms748tester<double>::reset();
 165 }
 166
 167 BOOST_AUTO_TEST_CASE( test_main )
 168 {
 169    std::cout << std::setprecision(18);
 170    run_test(3.14/2, 3.14, 1);
 171
 172    for(int i = 1; i <= 10; i += 1)
 173    {
 174       run_test(i*i + 1e-9, (i+1)*(i+1) - 1e-9, 2);
 175    }
 176
 177    run_test(-9.0, 31.0, 3, -40.0, -1.0);
 178    run_test(-9.0, 31.0, 3, -100.0, -2.0);
 179    run_test(-9.0, 31.0, 3, -200.0, -3.0);
 180
 181    for(int n = 4; n <= 12; n += 2)
 182    {
 183       run_test(0.0, 5.0, 4, 0.2, double(n));
 184    }
 185    for(int n = 4; n <= 12; n += 2)
 186    {
 187       run_test(0.0, 5.0, 4, 1.0, double(n));
 188    }
 189    for(int n = 8; n <= 14; n += 2)
 190    {
 191       run_test(-0.95, 4.05, 4, 1.0, double(n));
 192    }
 193    run_test(0.0, 1.5, 5);
 194    for(int n = 1; n <= 5; ++n)
 195    {
 196       run_test(0.0, 1.0, 6, n);
 197    }
 198    for(int n = 20; n <= 100; n += 20)
 199    {
 200       run_test(0.0, 1.0, 6, n);
 201    }
 202    run_test(0.0, 1.0, 7, 5);
 203    run_test(0.0, 1.0, 7, 10);
 204    run_test(0.0, 1.0, 7, 20);
 205    run_test(0.0, 1.0, 8, 2);
 206    run_test(0.0, 1.0, 8, 5);
 207    run_test(0.0, 1.0, 8, 10);
 208    run_test(0.0, 1.0, 8, 15);
 209    run_test(0.0, 1.0, 8, 20);
 210    run_test(0.0, 1.0, 9, 1);
 211    run_test(0.0, 1.0, 9, 2);
 212    run_test(0.0, 1.0, 9, 4);
 213    run_test(0.0, 1.0, 9, 5);
 214    run_test(0.0, 1.0, 9, 8);
 215    run_test(0.0, 1.0, 9, 15);
 216    run_test(0.0, 1.0, 9, 20);
 217    run_test(0.0, 1.0, 10, 1);
 218    run_test(0.0, 1.0, 10, 5);
 219    run_test(0.0, 1.0, 10, 10);
 220    run_test(0.0, 1.0, 10, 15);
 221    run_test(0.0, 1.0, 10, 20);
 222
 223    run_test(0.01, 1.0, 11, 2);
 224    run_test(0.01, 1.0, 11, 5);
 225    run_test(0.01, 1.0, 11, 15);
 226    run_test(0.01, 1.0, 11, 20);
 227
 228    for(int n = 2; n <= 6; ++n)
 229       run_test(1.0, 100.0, 12, n);
 230    for(int n = 7; n <= 33; n+=2)
 231       run_test(1.0, 100.0, 12, n);
 232
 233    run_test(-1.0, 4.0, 13);
 234
 235    for(int n = 1; n <= 40; ++n)
 236       run_test(-1e4, 3.14/2, 14, n);
 237
 238    for(int n = 20; n <= 40; ++n)
 239       run_test(-1e4, 1e-4, 15, n);
 240
 241    for(int n = 100; n <= 1000; n+=100)
 242       run_test(-1e4, 1e-4, 15, n);
 243
 244    std::cout << "Total iterations consumed = " << total << std::endl;
 245    std::cout << "Total function invocations consumed = " << invocations << std::endl << std::endl;
 246
 247    BOOST_CHECK(invocations < 3150);
 248
 249    std::cout << std::setprecision(18);
 250
 251    for(int n = 5; n <= 100; n += 10)
 252       run_test(sqrt(double(n)), double(n+1), 16, (double)n, 0.4);
 253
 254    for(int n = 5; n <= 100; n += 10)
 255       run_test(double(n / 2), double(2*n), 17, (double)n, 0.4);
 256
 257
 258    for(int n = 4; n <= 12; n += 2)
 259    {
 260       boost::uintmax_t c = 1000;
 261       std::pair<double, double> r = bracket_and_solve_root(toms748tester<double>(4, 0.2, double(n)),
 262          2.0,
 263          2.0,
 264          true,
 265          boost::math::tools::eps_tolerance<double>(std::numeric_limits<double>::digits),
 266          c);
 267       std::cout << std::setprecision(18);
 268       std::cout << "Function " << 4 << "\n   Result={" << r.first << ", " << r.second << "} total calls=" << toms748tester<double>::total_calls() << "\n\n";
 269       toms748tester<double>::reset();
 270       BOOST_CHECK(c < 20);
 271    }
 272 }
 273