]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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> | |
92f5a8d4 | 10 | #include <boost/test/tools/floating_point_comparison.hpp> |
7c673cae FG |
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 |