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 | // |
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(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 | } |
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(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 | |