]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // (C) Copyright John Maddock 2015. |
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 | #ifndef BOOST_NO_CXX11_HDR_TUPLE | |
9 | ||
10 | #define BOOST_TEST_MAIN | |
11 | #include <boost/test/unit_test.hpp> | |
12 | #include <boost/test/floating_point_comparison.hpp> | |
13 | #include <boost/math/tools/roots.hpp> | |
14 | #include <boost/test/results_collector.hpp> | |
15 | #include <boost/test/unit_test.hpp> | |
16 | #include <boost/math/special_functions/cbrt.hpp> | |
17 | #include <iostream> | |
18 | #include <iomanip> | |
19 | #include <tuple> | |
20 | ||
21 | // No derivatives - using TOMS748 internally. | |
22 | struct cbrt_functor_noderiv | |
23 | { // cube root of x using only function - no derivatives. | |
24 | cbrt_functor_noderiv(double to_find_root_of) : a(to_find_root_of) | |
25 | { // Constructor just stores value a to find root of. | |
26 | } | |
27 | double operator()(double x) | |
28 | { | |
29 | double fx = x*x*x - a; // Difference (estimate x^3 - a). | |
30 | return fx; | |
31 | } | |
32 | private: | |
33 | double a; // to be 'cube_rooted'. | |
34 | }; // template <class T> struct cbrt_functor_noderiv | |
35 | ||
36 | // Using 1st derivative only Newton-Raphson | |
37 | struct cbrt_functor_deriv | |
38 | { // Functor also returning 1st derviative. | |
39 | cbrt_functor_deriv(double const& to_find_root_of) : a(to_find_root_of) | |
40 | { // Constructor stores value a to find root of, | |
41 | // for example: calling cbrt_functor_deriv<double>(x) to use to get cube root of x. | |
42 | } | |
43 | std::pair<double, double> operator()(double const& x) | |
44 | { // Return both f(x) and f'(x). | |
45 | double fx = x*x*x - a; // Difference (estimate x^3 - value). | |
46 | double dx = 3 * x*x; // 1st derivative = 3x^2. | |
47 | return std::make_pair(fx, dx); // 'return' both fx and dx. | |
48 | } | |
49 | private: | |
50 | double a; // to be 'cube_rooted'. | |
51 | }; | |
52 | // Using 1st and 2nd derivatives with Halley algorithm. | |
53 | struct cbrt_functor_2deriv | |
54 | { // Functor returning both 1st and 2nd derivatives. | |
55 | cbrt_functor_2deriv(double const& to_find_root_of) : a(to_find_root_of) | |
56 | { // Constructor stores value a to find root of, for example: | |
57 | // calling cbrt_functor_2deriv<double>(x) to get cube root of x, | |
58 | } | |
59 | std::tuple<double, double, double> operator()(double const& x) | |
60 | { // Return both f(x) and f'(x) and f''(x). | |
61 | double fx = x*x*x - a; // Difference (estimate x^3 - value). | |
62 | double dx = 3 * x*x; // 1st derivative = 3x^2. | |
63 | double d2x = 6 * x; // 2nd derivative = 6x. | |
64 | return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x. | |
65 | } | |
66 | private: | |
67 | double a; // to be 'cube_rooted'. | |
68 | }; | |
69 | ||
70 | BOOST_AUTO_TEST_CASE( test_main ) | |
71 | { | |
72 | int newton_limits = static_cast<int>(std::numeric_limits<double>::digits * 0.6); | |
73 | int halley_limits = static_cast<int>(std::numeric_limits<double>::digits * 0.4); | |
74 | double arg = 1e-50; | |
75 | while(arg < 1e50) | |
76 | { | |
77 | double result = boost::math::cbrt(arg); | |
78 | // | |
79 | // Start with a really bad guess 5 times below the result: | |
80 | // | |
81 | double guess = result / 5; | |
82 | boost::uintmax_t iters = 1000; | |
83 | // TOMS algo first: | |
84 | std::pair<double, double> r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance<double>(), iters); | |
85 | BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4); | |
86 | BOOST_CHECK_LE(iters, 14); | |
87 | // Newton next: | |
88 | iters = 1000; | |
89 | double dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, guess / 2, result * 10, newton_limits, iters); | |
90 | BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2); | |
91 | BOOST_CHECK_LE(iters, 12); | |
92 | // Halley next: | |
93 | iters = 1000; | |
94 | dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters); | |
95 | BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2); | |
96 | BOOST_CHECK_LE(iters, 7); | |
97 | // Schroder next: | |
98 | iters = 1000; | |
99 | dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters); | |
100 | BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2); | |
101 | BOOST_CHECK_LE(iters, 11); | |
102 | // | |
103 | // Over again with a bad guess 5 times larger than the result: | |
104 | // | |
105 | iters = 1000; | |
106 | guess = result * 5; | |
107 | r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance<double>(), iters); | |
108 | BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4); | |
109 | BOOST_CHECK_LE(iters, 14); | |
110 | // Newton next: | |
111 | iters = 1000; | |
112 | dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters); | |
113 | BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2); | |
114 | BOOST_CHECK_LE(iters, 12); | |
115 | // Halley next: | |
116 | iters = 1000; | |
117 | dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters); | |
118 | BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2); | |
119 | BOOST_CHECK_LE(iters, 7); | |
120 | // Schroder next: | |
121 | iters = 1000; | |
122 | dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters); | |
123 | BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2); | |
124 | BOOST_CHECK_LE(iters, 11); | |
125 | // | |
126 | // A much better guess, 1% below result: | |
127 | // | |
128 | iters = 1000; | |
129 | guess = result * 0.9; | |
130 | r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance<double>(), iters); | |
131 | BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4); | |
132 | BOOST_CHECK_LE(iters, 12); | |
133 | // Newton next: | |
134 | iters = 1000; | |
135 | dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters); | |
136 | BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2); | |
137 | BOOST_CHECK_LE(iters, 5); | |
138 | // Halley next: | |
139 | iters = 1000; | |
140 | dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters); | |
141 | BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2); | |
142 | BOOST_CHECK_LE(iters, 3); | |
143 | // Schroder next: | |
144 | iters = 1000; | |
145 | dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters); | |
146 | BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2); | |
147 | BOOST_CHECK_LE(iters, 4); | |
148 | // | |
149 | // A much better guess, 1% above result: | |
150 | // | |
151 | iters = 1000; | |
152 | guess = result * 1.1; | |
153 | r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance<double>(), iters); | |
154 | BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits<double>::epsilon() * 4); | |
155 | BOOST_CHECK_LE(iters, 12); | |
156 | // Newton next: | |
157 | iters = 1000; | |
158 | dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters); | |
159 | BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2); | |
160 | BOOST_CHECK_LE(iters, 5); | |
161 | // Halley next: | |
162 | iters = 1000; | |
163 | dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters); | |
164 | BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2); | |
165 | BOOST_CHECK_LE(iters, 3); | |
166 | // Schroder next: | |
167 | iters = 1000; | |
168 | dr = boost::math::tools::schroder_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters); | |
169 | BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits<double>::epsilon() * 2); | |
170 | BOOST_CHECK_LE(iters, 4); | |
171 | ||
172 | arg *= 3.5; | |
173 | } | |
174 | } | |
175 | ||
176 | #else | |
177 | ||
178 | int main() { return 0; } | |
179 | ||
180 | #endif |