1 // Copyright Paul A. Bristow 2015.
2
3 // Use, modification and distribution are subject to the
4 // Boost Software License, Version 1.0.
5 // (See accompanying file LICENSE_1_0.txt
6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
7
8 // Note that this file contains Quickbook mark-up as well as code
9 // and comments, don't change any of the special comment mark-ups!
10
11 // Example of root finding using Boost.Multiprecision.
12
13 #include <boost/math/tools/roots.hpp>
14 //using boost::math::policies::policy;
15 //using boost::math::tools::newton_raphson_iterate;
16 //using boost::math::tools::halley_iterate;
17 //using boost::math::tools::eps_tolerance; // Binary functor for specified number of bits.
18 //using boost::math::tools::bracket_and_solve_root;
19 //using boost::math::tools::toms748_solve;
20
21 #include <boost/math/special_functions/next.hpp> // For float_distance.
22 #include <boost/math/special_functions/pow.hpp>
23 #include <boost/math/constants/constants.hpp>
24
25 //[root_finding_multiprecision_include_1
26 #include <boost/multiprecision/cpp_bin_float.hpp> // For cpp_bin_float_50.
27 #include <boost/multiprecision/cpp_dec_float.hpp> // For cpp_dec_float_50.
28 #ifndef _MSC_VER // float128 is not yet supported by Microsoft compiler at 2013.
29 # include <boost/multiprecision/float128.hpp> // Requires libquadmath.
30 #endif
31 //] [/root_finding_multiprecision_include_1]
32
33 #include <iostream>
34 // using std::cout; using std::endl;
35 #include <iomanip>
36 // using std::setw; using std::setprecision;
37 #include <limits>
38 // using std::numeric_limits;
39 #include <tuple>
40 #include <utility> // pair, make_pair
41
42 // #define BUILTIN_POW_GUESS // define to use std::pow function to obtain a guess.
43
44 template <class T>
cbrt_2deriv(T x)45 T cbrt_2deriv(T x)
46 { // return cube root of x using 1st and 2nd derivatives and Halley.
47 using namespace std; // Help ADL of std functions.
48 using namespace boost::math::tools; // For halley_iterate.
49
50 // If T is not a binary floating-point type, for example, cpp_dec_float_50
51 // then frexp may not be defined,
52 // so it may be necessary to compute the guess using a built-in type,
53 // probably quickest using double, but perhaps with float or long double.
54 // Note that the range of exponent may be restricted by a built-in-type for guess.
55
56 typedef long double guess_type;
57
58 #ifdef BUILTIN_POW_GUESS
59 guess_type pow_guess = std::pow(static_cast<guess_type>(x), static_cast<guess_type>(1) / 3);
60 T guess = pow_guess;
61 T min = pow_guess /2;
62 T max = pow_guess * 2;
63 #else
64 int exponent;
65 frexp(static_cast<guess_type>(x), &exponent); // Get exponent of z (ignore mantissa).
66 T guess = ldexp(static_cast<guess_type>(1.), exponent / 3); // Rough guess is to divide the exponent by three.
67 T min = ldexp(static_cast<guess_type>(1.) / 2, exponent / 3); // Minimum possible value is half our guess.
68 T max = ldexp(static_cast<guess_type>(2.), exponent / 3); // Maximum possible value is twice our guess.
69 #endif
70
71 int digits = std::numeric_limits<T>::digits / 2; // Half maximum possible binary digits accuracy for type T.
72 const boost::uintmax_t maxit = 20;
73 boost::uintmax_t it = maxit;
74 T result = halley_iterate(cbrt_functor_2deriv<T>(x), guess, min, max, digits, it);
75 // Can show how many iterations (updated by halley_iterate).
76 // std::cout << "Iterations " << it << " (from max of "<< maxit << ")." << std::endl;
77 return result;
78 } // cbrt_2deriv(x)
79
80
81 template <class T>
82 struct cbrt_functor_2deriv
83 { // Functor returning both 1st and 2nd derivatives.
cbrt_functor_2derivcbrt_functor_2deriv84 cbrt_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of)
85 { // Constructor stores value to find root of, for example:
86 }
87
88 // using boost::math::tuple; // to return three values.
operator ()cbrt_functor_2deriv89 std::tuple<T, T, T> operator()(T const& x)
90 {
91 // Return both f(x) and f'(x) and f''(x).
92 T fx = x*x*x - a; // Difference (estimate x^3 - value).
93 // std::cout << "x = " << x << "\nfx = " << fx << std::endl;
94 T dx = 3 * x*x; // 1st derivative = 3x^2.
95 T d2x = 6 * x; // 2nd derivative = 6x.
96 return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x.
97 }
98 private:
99 T a; // to be 'cube_rooted'.
100 }; // struct cbrt_functor_2deriv
101
102 template <int n, class T>
103 struct nth_functor_2deriv
104 { // Functor returning both 1st and 2nd derivatives.
105
nth_functor_2derivnth_functor_2deriv106 nth_functor_2deriv(T const& to_find_root_of) : value(to_find_root_of)
107 { /* Constructor stores value to find root of, for example: */ }
108
109 // using std::tuple; // to return three values.
operator ()nth_functor_2deriv110 std::tuple<T, T, T> operator()(T const& x)
111 {
112 // Return both f(x) and f'(x) and f''(x).
113 using boost::math::pow;
114 T fx = pow<n>(x) - value; // Difference (estimate x^3 - value).
115 T dx = n * pow<n - 1>(x); // 1st derivative = 5x^4.
116 T d2x = n * (n - 1) * pow<n - 2 >(x); // 2nd derivative = 20 x^3
117 return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x.
118 }
119 private:
120 T value; // to be 'nth_rooted'.
121 }; // struct nth_functor_2deriv
122
123
124 template <int n, class T>
nth_2deriv(T x)125 T nth_2deriv(T x)
126 {
127 // return nth root of x using 1st and 2nd derivatives and Halley.
128 using namespace std; // Help ADL of std functions.
129 using namespace boost::math; // For halley_iterate.
130
131 int exponent;
132 frexp(x, &exponent); // Get exponent of z (ignore mantissa).
133 T guess = ldexp(static_cast<T>(1.), exponent / n); // Rough guess is to divide the exponent by three.
134 T min = ldexp(static_cast<T>(0.5), exponent / n); // Minimum possible value is half our guess.
135 T max = ldexp(static_cast<T>(2.), exponent / n); // Maximum possible value is twice our guess.
136
137 int digits = std::numeric_limits<T>::digits / 2; // Half maximum possible binary digits accuracy for type T.
138 const boost::uintmax_t maxit = 50;
139 boost::uintmax_t it = maxit;
140 T result = halley_iterate(nth_functor_2deriv<n, T>(x), guess, min, max, digits, it);
141 // Can show how many iterations (updated by halley_iterate).
142 std::cout << it << " iterations (from max of " << maxit << ")" << std::endl;
143
144 return result;
145 } // nth_2deriv(x)
146
147 //[root_finding_multiprecision_show_1
148
149 template <typename T>
show_cube_root(T value)150 T show_cube_root(T value)
151 { // Demonstrate by printing the root using all definitely significant digits.
152 std::cout.precision(std::numeric_limits<T>::digits10);
153 T r = cbrt_2deriv(value);
154 std::cout << "value = " << value << ", cube root =" << r << std::endl;
155 return r;
156 }
157
158 //] [/root_finding_multiprecision_show_1]
159
main()160 int main()
161 {
162 std::cout << "Multiprecision Root finding Example." << std::endl;
163 // Show all possibly significant decimal digits.
164 std::cout.precision(std::numeric_limits<double>::digits10);
165 // or use cout.precision(max_digits10 = 2 + std::numeric_limits<double>::digits * 3010/10000);
166 //[root_finding_multiprecision_example_1
167 using boost::multiprecision::cpp_dec_float_50; // decimal.
168 using boost::multiprecision::cpp_bin_float_50; // binary.
169 #ifndef _MSC_VER // Not supported by Microsoft compiler.
170 using boost::multiprecision::float128;
171 #endif
172 //] [/root_finding_multiprecision_example_1
173
174 try
175 { // Always use try'n'catch blocks with Boost.Math to get any error messages.
176 // Increase the precision to 50 decimal digits using Boost.Multiprecision
177 //[root_finding_multiprecision_example_2
178
179 std::cout.precision(std::numeric_limits<cpp_dec_float_50>::digits10);
180
181 cpp_dec_float_50 two = 2; //
182 cpp_dec_float_50 r = cbrt_2deriv(two);
183 std::cout << "cbrt(" << two << ") = " << r << std::endl;
184
185 r = cbrt_2deriv(2.); // Passing a double, so ADL will compute a double precision result.
186 std::cout << "cbrt(" << two << ") = " << r << std::endl;
187 // cbrt(2) = 1.2599210498948731906665443602832965552806854248047 'wrong' from digits 17 onwards!
188 r = cbrt_2deriv(static_cast<cpp_dec_float_50>(2.)); // Passing a cpp_dec_float_50,
189 // so will compute a cpp_dec_float_50 precision result.
190 std::cout << "cbrt(" << two << ") = " << r << std::endl;
191 r = cbrt_2deriv<cpp_dec_float_50>(2.); // Explictly a cpp_dec_float_50, so will compute a cpp_dec_float_50 precision result.
192 std::cout << "cbrt(" << two << ") = " << r << std::endl;
193 // cpp_dec_float_50 1.2599210498948731647672106072782283505702514647015
194 //] [/root_finding_multiprecision_example_2
195 // N[2^(1/3), 50] 1.2599210498948731647672106072782283505702514647015
196
197 //show_cube_root(2); // Integer parameter - Errors!
198 //show_cube_root(2.F); // Float parameter - Warnings!
199 //[root_finding_multiprecision_example_3
200 show_cube_root(2.);
201 show_cube_root(2.L);
202 show_cube_root(two);
203
204 //] [/root_finding_multiprecision_example_3
205
206 }
207 catch (const std::exception& e)
208 { // Always useful to include try&catch blocks because default policies
209 // are to throw exceptions on arguments that cause errors like underflow & overflow.
210 // Lacking try&catch blocks, the program will abort without a message below,
211 // which may give some helpful clues as to the cause of the exception.
212 std::cout <<
213 "\n""Message from thrown exception was:\n " << e.what() << std::endl;
214 }
215 return 0;
216 } // int main()
217
218
219 /*
220
221 Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Release\root_finding_multiprecision.exe"
222 Multiprecision Root finding Example.
223 cbrt(2) = 1.2599210498948731647672106072782283505702514647015
224 cbrt(2) = 1.2599210498948731906665443602832965552806854248047
225 cbrt(2) = 1.2599210498948731647672106072782283505702514647015
226 cbrt(2) = 1.2599210498948731647672106072782283505702514647015
227 value = 2, cube root =1.25992104989487
228 value = 2, cube root =1.25992104989487
229 value = 2, cube root =1.2599210498948731647672106072782283505702514647015
230
231
232 */
233