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/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//
26template <class T>
27struct 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(x: 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 }
105private:
106 int k; // index of problem.
107 int ip; // integer parameter
108 T a, b; // real parameter
109
110 static unsigned invocations;
111};
112
113template <class T>
114unsigned toms748tester<T>::invocations = 0;
115
116boost::uintmax_t total = 0;
117boost::uintmax_t invocations = 0;
118
119template <class T>
120void 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
135template <class T>
136void 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
151template <class T>
152void 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
167BOOST_AUTO_TEST_CASE( test_main )
168{
169 std::cout << std::setprecision(18);
170 run_test(a: 3.14/2, b: 3.14, id: 1);
171
172 for(int i = 1; i <= 10; i += 1)
173 {
174 run_test(a: i*i + 1e-9, b: (i+1)*(i+1) - 1e-9, id: 2);
175 }
176
177 run_test(a: -9.0, b: 31.0, id: 3, p1: -40.0, p2: -1.0);
178 run_test(a: -9.0, b: 31.0, id: 3, p1: -100.0, p2: -2.0);
179 run_test(a: -9.0, b: 31.0, id: 3, p1: -200.0, p2: -3.0);
180
181 for(int n = 4; n <= 12; n += 2)
182 {
183 run_test(a: 0.0, b: 5.0, id: 4, p1: 0.2, p2: double(n));
184 }
185 for(int n = 4; n <= 12; n += 2)
186 {
187 run_test(a: 0.0, b: 5.0, id: 4, p1: 1.0, p2: double(n));
188 }
189 for(int n = 8; n <= 14; n += 2)
190 {
191 run_test(a: -0.95, b: 4.05, id: 4, p1: 1.0, p2: double(n));
192 }
193 run_test(a: 0.0, b: 1.5, id: 5);
194 for(int n = 1; n <= 5; ++n)
195 {
196 run_test(a: 0.0, b: 1.0, id: 6, p: n);
197 }
198 for(int n = 20; n <= 100; n += 20)
199 {
200 run_test(a: 0.0, b: 1.0, id: 6, p: n);
201 }
202 run_test(a: 0.0, b: 1.0, id: 7, p: 5);
203 run_test(a: 0.0, b: 1.0, id: 7, p: 10);
204 run_test(a: 0.0, b: 1.0, id: 7, p: 20);
205 run_test(a: 0.0, b: 1.0, id: 8, p: 2);
206 run_test(a: 0.0, b: 1.0, id: 8, p: 5);
207 run_test(a: 0.0, b: 1.0, id: 8, p: 10);
208 run_test(a: 0.0, b: 1.0, id: 8, p: 15);
209 run_test(a: 0.0, b: 1.0, id: 8, p: 20);
210 run_test(a: 0.0, b: 1.0, id: 9, p: 1);
211 run_test(a: 0.0, b: 1.0, id: 9, p: 2);
212 run_test(a: 0.0, b: 1.0, id: 9, p: 4);
213 run_test(a: 0.0, b: 1.0, id: 9, p: 5);
214 run_test(a: 0.0, b: 1.0, id: 9, p: 8);
215 run_test(a: 0.0, b: 1.0, id: 9, p: 15);
216 run_test(a: 0.0, b: 1.0, id: 9, p: 20);
217 run_test(a: 0.0, b: 1.0, id: 10, p: 1);
218 run_test(a: 0.0, b: 1.0, id: 10, p: 5);
219 run_test(a: 0.0, b: 1.0, id: 10, p: 10);
220 run_test(a: 0.0, b: 1.0, id: 10, p: 15);
221 run_test(a: 0.0, b: 1.0, id: 10, p: 20);
222
223 run_test(a: 0.01, b: 1.0, id: 11, p: 2);
224 run_test(a: 0.01, b: 1.0, id: 11, p: 5);
225 run_test(a: 0.01, b: 1.0, id: 11, p: 15);
226 run_test(a: 0.01, b: 1.0, id: 11, p: 20);
227
228 for(int n = 2; n <= 6; ++n)
229 run_test(a: 1.0, b: 100.0, id: 12, p: n);
230 for(int n = 7; n <= 33; n+=2)
231 run_test(a: 1.0, b: 100.0, id: 12, p: n);
232
233 run_test(a: -1.0, b: 4.0, id: 13);
234
235 for(int n = 1; n <= 40; ++n)
236 run_test(a: -1e4, b: 3.14/2, id: 14, p: n);
237
238 for(int n = 20; n <= 40; ++n)
239 run_test(a: -1e4, b: 1e-4, id: 15, p: n);
240
241 for(int n = 100; n <= 1000; n+=100)
242 run_test(a: -1e4, b: 1e-4, id: 15, p: 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(a: sqrt(x: double(n)), b: double(n+1), id: 16, p1: (double)n, p2: 0.4);
253
254 for(int n = 5; n <= 100; n += 10)
255 run_test(a: double(n / 2), b: double(2*n), id: 17, p1: (double)n, p2: 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(f: toms748tester<double>(4, 0.2, double(n)),
262 guess: 2.0,
263 factor: 2.0,
264 rising: true,
265 tol: boost::math::tools::eps_tolerance<double>(std::numeric_limits<double>::digits),
266 max_iter&: 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
274

source code of boost/libs/math/test/test_toms748_solve.cpp