1 // Copyright John Maddock 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 // Comparison of finding roots using TOMS748, Newton-Raphson, Halley & Schroder algorithms.
9 // Note that this file contains Quickbook mark-up as well as code
10 // and comments, don't change any of the special comment mark-ups!
11 // This program also writes files in Quickbook tables mark-up format.
12
13 #include <boost/cstdlib.hpp>
14 #include <boost/config.hpp>
15 #include <boost/array.hpp>
16 #include <boost/math/tools/roots.hpp>
17 #include <boost/math/special_functions/ellint_1.hpp>
18 #include <boost/math/special_functions/ellint_2.hpp>
19 template <class T>
20 struct cbrt_functor_noderiv
21 {
22 // cube root of x using only function - no derivatives.
cbrt_functor_noderivcbrt_functor_noderiv23 cbrt_functor_noderiv(T const& to_find_root_of) : a(to_find_root_of)
24 { /* Constructor just stores value a to find root of. */
25 }
operator ()cbrt_functor_noderiv26 T operator()(T const& x)
27 {
28 T fx = x*x*x - a; // Difference (estimate x^3 - a).
29 return fx;
30 }
31 private:
32 T a; // to be 'cube_rooted'.
33 };
34 //] [/root_finding_noderiv_1
35
36 template <class T>
cbrt_noderiv(T x,T guess)37 boost::uintmax_t cbrt_noderiv(T x, T guess)
38 {
39 // return cube root of x using bracket_and_solve (no derivatives).
40 using namespace std; // Help ADL of std functions.
41 using namespace boost::math::tools; // For bracket_and_solve_root.
42
43 T factor = 2; // How big steps to take when searching.
44
45 const boost::uintmax_t maxit = 20; // Limit to maximum iterations.
46 boost::uintmax_t it = maxit; // Initally our chosen max iterations, but updated with actual.
47 bool is_rising = true; // So if result if guess^3 is too low, then try increasing guess.
48 int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.
49 // Some fraction of digits is used to control how accurate to try to make the result.
50 int get_digits = digits - 3; // We have to have a non-zero interval at each step, so
51 // maximum accuracy is digits - 1. But we also have to
52 // allow for inaccuracy in f(x), otherwise the last few
53 // iterations just thrash around.
54 eps_tolerance<T> tol(get_digits); // Set the tolerance.
55 bracket_and_solve_root(cbrt_functor_noderiv<T>(x), guess, factor, is_rising, tol, it);
56 return it;
57 }
58
59 template <class T>
60 struct cbrt_functor_deriv
61 { // Functor also returning 1st derivative.
cbrt_functor_derivcbrt_functor_deriv62 cbrt_functor_deriv(T const& to_find_root_of) : a(to_find_root_of)
63 { // Constructor stores value a to find root of,
64 // for example: calling cbrt_functor_deriv<T>(a) to use to get cube root of a.
65 }
operator ()cbrt_functor_deriv66 std::pair<T, T> operator()(T const& x)
67 {
68 // Return both f(x) and f'(x).
69 T fx = x*x*x - a; // Difference (estimate x^3 - value).
70 T dx = 3 * x*x; // 1st derivative = 3x^2.
71 return std::make_pair(fx, dx); // 'return' both fx and dx.
72 }
73 private:
74 T a; // Store value to be 'cube_rooted'.
75 };
76
77 template <class T>
cbrt_deriv(T x,T guess)78 boost::uintmax_t cbrt_deriv(T x, T guess)
79 {
80 // return cube root of x using 1st derivative and Newton_Raphson.
81 using namespace boost::math::tools;
82 T min = guess / 100; // We don't really know what this should be!
83 T max = guess * 100; // We don't really know what this should be!
84 const int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.
85 int get_digits = static_cast<int>(digits * 0.6); // Accuracy doubles with each step, so stop when we have
86 // just over half the digits correct.
87 const boost::uintmax_t maxit = 20;
88 boost::uintmax_t it = maxit;
89 newton_raphson_iterate(cbrt_functor_deriv<T>(x), guess, min, max, get_digits, it);
90 return it;
91 }
92
93 template <class T>
94 struct cbrt_functor_2deriv
95 {
96 // Functor returning both 1st and 2nd derivatives.
cbrt_functor_2derivcbrt_functor_2deriv97 cbrt_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of)
98 { // Constructor stores value a to find root of, for example:
99 // calling cbrt_functor_2deriv<T>(x) to get cube root of x,
100 }
operator ()cbrt_functor_2deriv101 std::tuple<T, T, T> operator()(T const& x)
102 {
103 // Return both f(x) and f'(x) and f''(x).
104 T fx = x*x*x - a; // Difference (estimate x^3 - value).
105 T dx = 3 * x*x; // 1st derivative = 3x^2.
106 T d2x = 6 * x; // 2nd derivative = 6x.
107 return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x.
108 }
109 private:
110 T a; // to be 'cube_rooted'.
111 };
112
113 template <class T>
cbrt_2deriv(T x,T guess)114 boost::uintmax_t cbrt_2deriv(T x, T guess)
115 {
116 // return cube root of x using 1st and 2nd derivatives and Halley.
117 //using namespace std; // Help ADL of std functions.
118 using namespace boost::math::tools;
119 T min = guess / 100; // We don't really know what this should be!
120 T max = guess * 100; // We don't really know what this should be!
121 const int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.
122 // digits used to control how accurate to try to make the result.
123 int get_digits = static_cast<int>(digits * 0.4); // Accuracy triples with each step, so stop when just
124 // over one third of the digits are correct.
125 boost::uintmax_t maxit = 20;
126 halley_iterate(cbrt_functor_2deriv<T>(x), guess, min, max, get_digits, maxit);
127 return maxit;
128 }
129
130 template <class T>
cbrt_2deriv_s(T x,T guess)131 boost::uintmax_t cbrt_2deriv_s(T x, T guess)
132 {
133 // return cube root of x using 1st and 2nd derivatives and Halley.
134 //using namespace std; // Help ADL of std functions.
135 using namespace boost::math::tools;
136 T min = guess / 100; // We don't really know what this should be!
137 T max = guess * 100; // We don't really know what this should be!
138 const int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.
139 // digits used to control how accurate to try to make the result.
140 int get_digits = static_cast<int>(digits * 0.4); // Accuracy triples with each step, so stop when just
141 // over one third of the digits are correct.
142 boost::uintmax_t maxit = 20;
143 schroder_iterate(cbrt_functor_2deriv<T>(x), guess, min, max, get_digits, maxit);
144 return maxit;
145 }
146
147 template <typename T = double>
148 struct elliptic_root_functor_noderiv
149 {
elliptic_root_functor_noderivelliptic_root_functor_noderiv150 elliptic_root_functor_noderiv(T const& arc, T const& radius) : m_arc(arc), m_radius(radius)
151 { // Constructor just stores value a to find root of.
152 }
operator ()elliptic_root_functor_noderiv153 T operator()(T const& x)
154 {
155 // return the difference between required arc-length, and the calculated arc-length for an
156 // ellipse with radii m_radius and x:
157 T a = (std::max)(m_radius, x);
158 T b = (std::min)(m_radius, x);
159 T k = sqrt(1 - b * b / (a * a));
160 return 4 * a * boost::math::ellint_2(k) - m_arc;
161 }
162 private:
163 T m_arc; // length of arc.
164 T m_radius; // one of the two radii of the ellipse
165 }; // template <class T> struct elliptic_root_functor_noderiv
166
167 template <class T = double>
elliptic_root_noderiv(T radius,T arc,T guess)168 boost::uintmax_t elliptic_root_noderiv(T radius, T arc, T guess)
169 { // return the other radius of an ellipse, given one radii and the arc-length
170 using namespace std; // Help ADL of std functions.
171 using namespace boost::math::tools; // For bracket_and_solve_root.
172
173 T factor = 2; // How big steps to take when searching.
174
175 const boost::uintmax_t maxit = 50; // Limit to maximum iterations.
176 boost::uintmax_t it = maxit; // Initally our chosen max iterations, but updated with actual.
177 bool is_rising = true; // arc-length increases if one radii increases, so function is rising
178 // Define a termination condition, stop when nearly all digits are correct, but allow for
179 // the fact that we are returning a range, and must have some inaccuracy in the elliptic integral:
180 eps_tolerance<T> tol(std::numeric_limits<T>::digits - 2);
181 // Call bracket_and_solve_root to find the solution, note that this is a rising function:
182 bracket_and_solve_root(elliptic_root_functor_noderiv<T>(arc, radius), guess, factor, is_rising, tol, it);
183 return it;
184 }
185
186 template <class T = double>
187 struct elliptic_root_functor_1deriv
188 { // Functor also returning 1st derviative.
189 BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!");
190
elliptic_root_functor_1derivelliptic_root_functor_1deriv191 elliptic_root_functor_1deriv(T const& arc, T const& radius) : m_arc(arc), m_radius(radius)
192 { // Constructor just stores value a to find root of.
193 }
operator ()elliptic_root_functor_1deriv194 std::pair<T, T> operator()(T const& x)
195 {
196 // Return the difference between required arc-length, and the calculated arc-length for an
197 // ellipse with radii m_radius and x, plus it's derivative.
198 // See http://www.wolframalpha.com/input/?i=d%2Fda+[4+*+a+*+EllipticE%281+-+b^2%2Fa^2%29]
199 // We require two elliptic integral calls, but from these we can calculate both
200 // the function and it's derivative:
201 T a = (std::max)(m_radius, x);
202 T b = (std::min)(m_radius, x);
203 T a2 = a * a;
204 T b2 = b * b;
205 T k = sqrt(1 - b2 / a2);
206 T Ek = boost::math::ellint_2(k);
207 T Kk = boost::math::ellint_1(k);
208 T fx = 4 * a * Ek - m_arc;
209 T dfx = 4 * (a2 * Ek - b2 * Kk) / (a2 - b2);
210 return std::make_pair(fx, dfx);
211 }
212 private:
213 T m_arc; // length of arc.
214 T m_radius; // one of the two radii of the ellipse
215 }; // struct elliptic_root__functor_1deriv
216
217 template <class T = double>
elliptic_root_1deriv(T radius,T arc,T guess)218 boost::uintmax_t elliptic_root_1deriv(T radius, T arc, T guess)
219 {
220 using namespace std; // Help ADL of std functions.
221 using namespace boost::math::tools; // For newton_raphson_iterate.
222
223 BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!");
224
225 T min = 0; // Minimum possible value is zero.
226 T max = arc; // Maximum possible value is the arc length.
227
228 // Accuracy doubles at each step, so stop when just over half of the digits are
229 // correct, and rely on that step to polish off the remainder:
230 int get_digits = static_cast<int>(std::numeric_limits<T>::digits * 0.6);
231 const boost::uintmax_t maxit = 20;
232 boost::uintmax_t it = maxit;
233 newton_raphson_iterate(elliptic_root_functor_1deriv<T>(arc, radius), guess, min, max, get_digits, it);
234 return it;
235 }
236
237 template <class T = double>
238 struct elliptic_root_functor_2deriv
239 { // Functor returning both 1st and 2nd derivatives.
240 BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!");
241
elliptic_root_functor_2derivelliptic_root_functor_2deriv242 elliptic_root_functor_2deriv(T const& arc, T const& radius) : m_arc(arc), m_radius(radius) {}
operator ()elliptic_root_functor_2deriv243 std::tuple<T, T, T> operator()(T const& x)
244 {
245 // Return the difference between required arc-length, and the calculated arc-length for an
246 // ellipse with radii m_radius and x, plus it's derivative.
247 // See http://www.wolframalpha.com/input/?i=d^2%2Fda^2+[4+*+a+*+EllipticE%281+-+b^2%2Fa^2%29]
248 // for the second derivative.
249 T a = (std::max)(m_radius, x);
250 T b = (std::min)(m_radius, x);
251 T a2 = a * a;
252 T b2 = b * b;
253 T k = sqrt(1 - b2 / a2);
254 T Ek = boost::math::ellint_2(k);
255 T Kk = boost::math::ellint_1(k);
256 T fx = 4 * a * Ek - m_arc;
257 T dfx = 4 * (a2 * Ek - b2 * Kk) / (a2 - b2);
258 T dfx2 = 4 * b2 * ((a2 + b2) * Kk - 2 * a2 * Ek) / (a * (a2 - b2) * (a2 - b2));
259 return std::make_tuple(fx, dfx, dfx2);
260 }
261 private:
262 T m_arc; // length of arc.
263 T m_radius; // one of the two radii of the ellipse
264 };
265
266 template <class T = double>
elliptic_root_2deriv(T radius,T arc,T guess)267 boost::uintmax_t elliptic_root_2deriv(T radius, T arc, T guess)
268 {
269 using namespace std; // Help ADL of std functions.
270 using namespace boost::math::tools; // For halley_iterate.
271
272 BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!");
273
274 T min = 0; // Minimum possible value is zero.
275 T max = arc; // radius can't be larger than the arc length.
276
277 // Accuracy triples at each step, so stop when just over one-third of the digits
278 // are correct, and the last iteration will polish off the remaining digits:
279 int get_digits = static_cast<int>(std::numeric_limits<T>::digits * 0.4);
280 const boost::uintmax_t maxit = 20;
281 boost::uintmax_t it = maxit;
282 halley_iterate(elliptic_root_functor_2deriv<T>(arc, radius), guess, min, max, get_digits, it);
283 return it;
284 } // nth_2deriv Halley
285 //]
286 // Using 1st and 2nd derivatives using Schroder algorithm.
287
288 template <class T = double>
elliptic_root_2deriv_s(T radius,T arc,T guess)289 boost::uintmax_t elliptic_root_2deriv_s(T radius, T arc, T guess)
290 { // return nth root of x using 1st and 2nd derivatives and Schroder.
291
292 using namespace std; // Help ADL of std functions.
293 using namespace boost::math::tools; // For schroder_iterate.
294
295 BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!");
296
297 T min = 0; // Minimum possible value is zero.
298 T max = arc; // radius can't be larger than the arc length.
299
300 int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.
301 int get_digits = static_cast<int>(digits * 0.4);
302 const boost::uintmax_t maxit = 20;
303 boost::uintmax_t it = maxit;
304 schroder_iterate(elliptic_root_functor_2deriv<T>(arc, radius), guess, min, max, get_digits, it);
305 return it;
306 } // T elliptic_root_2deriv_s Schroder
307
308
main()309 int main()
310 {
311 try
312 {
313 double to_root = 500;
314 double answer = 7.93700525984;
315
316 std::cout << "[table\n"
317 << "[[Initial Guess=][-500% ([approx]1.323)][-100% ([approx]3.97)][-50% ([approx]3.96)][-20% ([approx]6.35)][-10% ([approx]7.14)][-5% ([approx]7.54)]"
318 "[5% ([approx]8.33)][10% ([approx]8.73)][20% ([approx]9.52)][50% ([approx]11.91)][100% ([approx]15.87)][500 ([approx]47.6)]]\n";
319 std::cout << "[[bracket_and_solve_root]["
320 << cbrt_noderiv(to_root, answer / 6)
321 << "][" << cbrt_noderiv(to_root, answer / 2)
322 << "][" << cbrt_noderiv(to_root, answer - answer * 0.5)
323 << "][" << cbrt_noderiv(to_root, answer - answer * 0.2)
324 << "][" << cbrt_noderiv(to_root, answer - answer * 0.1)
325 << "][" << cbrt_noderiv(to_root, answer - answer * 0.05)
326 << "][" << cbrt_noderiv(to_root, answer + answer * 0.05)
327 << "][" << cbrt_noderiv(to_root, answer + answer * 0.1)
328 << "][" << cbrt_noderiv(to_root, answer + answer * 0.2)
329 << "][" << cbrt_noderiv(to_root, answer + answer * 0.5)
330 << "][" << cbrt_noderiv(to_root, answer + answer)
331 << "][" << cbrt_noderiv(to_root, answer + answer * 5) << "]]\n";
332
333 std::cout << "[[newton_iterate]["
334 << cbrt_deriv(to_root, answer / 6)
335 << "][" << cbrt_deriv(to_root, answer / 2)
336 << "][" << cbrt_deriv(to_root, answer - answer * 0.5)
337 << "][" << cbrt_deriv(to_root, answer - answer * 0.2)
338 << "][" << cbrt_deriv(to_root, answer - answer * 0.1)
339 << "][" << cbrt_deriv(to_root, answer - answer * 0.05)
340 << "][" << cbrt_deriv(to_root, answer + answer * 0.05)
341 << "][" << cbrt_deriv(to_root, answer + answer * 0.1)
342 << "][" << cbrt_deriv(to_root, answer + answer * 0.2)
343 << "][" << cbrt_deriv(to_root, answer + answer * 0.5)
344 << "][" << cbrt_deriv(to_root, answer + answer)
345 << "][" << cbrt_deriv(to_root, answer + answer * 5) << "]]\n";
346
347 std::cout << "[[halley_iterate]["
348 << cbrt_2deriv(to_root, answer / 6)
349 << "][" << cbrt_2deriv(to_root, answer / 2)
350 << "][" << cbrt_2deriv(to_root, answer - answer * 0.5)
351 << "][" << cbrt_2deriv(to_root, answer - answer * 0.2)
352 << "][" << cbrt_2deriv(to_root, answer - answer * 0.1)
353 << "][" << cbrt_2deriv(to_root, answer - answer * 0.05)
354 << "][" << cbrt_2deriv(to_root, answer + answer * 0.05)
355 << "][" << cbrt_2deriv(to_root, answer + answer * 0.1)
356 << "][" << cbrt_2deriv(to_root, answer + answer * 0.2)
357 << "][" << cbrt_2deriv(to_root, answer + answer * 0.5)
358 << "][" << cbrt_2deriv(to_root, answer + answer)
359 << "][" << cbrt_2deriv(to_root, answer + answer * 5) << "]]\n";
360
361 std::cout << "[[schr'''ö'''der_iterate]["
362 << cbrt_2deriv_s(to_root, answer / 6)
363 << "][" << cbrt_2deriv_s(to_root, answer / 2)
364 << "][" << cbrt_2deriv_s(to_root, answer - answer * 0.5)
365 << "][" << cbrt_2deriv_s(to_root, answer - answer * 0.2)
366 << "][" << cbrt_2deriv_s(to_root, answer - answer * 0.1)
367 << "][" << cbrt_2deriv_s(to_root, answer - answer * 0.05)
368 << "][" << cbrt_2deriv_s(to_root, answer + answer * 0.05)
369 << "][" << cbrt_2deriv_s(to_root, answer + answer * 0.1)
370 << "][" << cbrt_2deriv_s(to_root, answer + answer * 0.2)
371 << "][" << cbrt_2deriv_s(to_root, answer + answer * 0.5)
372 << "][" << cbrt_2deriv_s(to_root, answer + answer)
373 << "][" << cbrt_2deriv_s(to_root, answer + answer * 5) << "]]\n]\n\n";
374
375
376 double radius_a = 10;
377 double arc_length = 500;
378 double radius_b = 123.6216507967705;
379
380 std::cout << std::setprecision(4) << "[table\n"
381 << "[[Initial Guess=][-500% ([approx]" << radius_b / 6 << ")][-100% ([approx]" << radius_b / 2 << ")][-50% ([approx]"
382 << radius_b - radius_b * 0.5 << ")][-20% ([approx]" << radius_b - radius_b * 0.2 << ")][-10% ([approx]" << radius_b - radius_b * 0.1 << ")][-5% ([approx]" << radius_b - radius_b * 0.05 << ")]"
383 "[5% ([approx]" << radius_b + radius_b * 0.05 << ")][10% ([approx]" << radius_b + radius_b * 0.1 << ")][20% ([approx]" << radius_b + radius_b * 0.2 << ")][50% ([approx]" << radius_b + radius_b * 0.5
384 << ")][100% ([approx]" << radius_b + radius_b << ")][500 ([approx]" << radius_b + radius_b * 5 << ")]]\n";
385 std::cout << "[[bracket_and_solve_root]["
386 << elliptic_root_noderiv(radius_a, arc_length, radius_b / 6)
387 << "][" << elliptic_root_noderiv(radius_a, arc_length, radius_b / 2)
388 << "][" << elliptic_root_noderiv(radius_a, arc_length, radius_b - radius_b * 0.5)
389 << "][" << elliptic_root_noderiv(radius_a, arc_length, radius_b - radius_b * 0.2)
390 << "][" << elliptic_root_noderiv(radius_a, arc_length, radius_b - radius_b * 0.1)
391 << "][" << elliptic_root_noderiv(radius_a, arc_length, radius_b - radius_b * 0.05)
392 << "][" << elliptic_root_noderiv(radius_a, arc_length, radius_b + radius_b * 0.05)
393 << "][" << elliptic_root_noderiv(radius_a, arc_length, radius_b + radius_b * 0.1)
394 << "][" << elliptic_root_noderiv(radius_a, arc_length, radius_b + radius_b * 0.2)
395 << "][" << elliptic_root_noderiv(radius_a, arc_length, radius_b + radius_b * 0.5)
396 << "][" << elliptic_root_noderiv(radius_a, arc_length, radius_b + radius_b)
397 << "][" << elliptic_root_noderiv(radius_a, arc_length, radius_b + radius_b * 5) << "]]\n";
398
399 std::cout << "[[newton_iterate]["
400 << elliptic_root_1deriv(radius_a, arc_length, radius_b / 6)
401 << "][" << elliptic_root_1deriv(radius_a, arc_length, radius_b / 2)
402 << "][" << elliptic_root_1deriv(radius_a, arc_length, radius_b - radius_b * 0.5)
403 << "][" << elliptic_root_1deriv(radius_a, arc_length, radius_b - radius_b * 0.2)
404 << "][" << elliptic_root_1deriv(radius_a, arc_length, radius_b - radius_b * 0.1)
405 << "][" << elliptic_root_1deriv(radius_a, arc_length, radius_b - radius_b * 0.05)
406 << "][" << elliptic_root_1deriv(radius_a, arc_length, radius_b + radius_b * 0.05)
407 << "][" << elliptic_root_1deriv(radius_a, arc_length, radius_b + radius_b * 0.1)
408 << "][" << elliptic_root_1deriv(radius_a, arc_length, radius_b + radius_b * 0.2)
409 << "][" << elliptic_root_1deriv(radius_a, arc_length, radius_b + radius_b * 0.5)
410 << "][" << elliptic_root_1deriv(radius_a, arc_length, radius_b + radius_b)
411 << "][" << elliptic_root_1deriv(radius_a, arc_length, radius_b + radius_b * 5) << "]]\n";
412
413 std::cout << "[[halley_iterate]["
414 << elliptic_root_2deriv(radius_a, arc_length, radius_b / 6)
415 << "][" << elliptic_root_2deriv(radius_a, arc_length, radius_b / 2)
416 << "][" << elliptic_root_2deriv(radius_a, arc_length, radius_b - radius_b * 0.5)
417 << "][" << elliptic_root_2deriv(radius_a, arc_length, radius_b - radius_b * 0.2)
418 << "][" << elliptic_root_2deriv(radius_a, arc_length, radius_b - radius_b * 0.1)
419 << "][" << elliptic_root_2deriv(radius_a, arc_length, radius_b - radius_b * 0.05)
420 << "][" << elliptic_root_2deriv(radius_a, arc_length, radius_b + radius_b * 0.05)
421 << "][" << elliptic_root_2deriv(radius_a, arc_length, radius_b + radius_b * 0.1)
422 << "][" << elliptic_root_2deriv(radius_a, arc_length, radius_b + radius_b * 0.2)
423 << "][" << elliptic_root_2deriv(radius_a, arc_length, radius_b + radius_b * 0.5)
424 << "][" << elliptic_root_2deriv(radius_a, arc_length, radius_b + radius_b)
425 << "][" << elliptic_root_2deriv(radius_a, arc_length, radius_b + radius_b * 5) << "]]\n";
426
427 std::cout << "[[schr'''ö'''der_iterate]["
428 << elliptic_root_2deriv_s(radius_a, arc_length, radius_b / 6)
429 << "][" << elliptic_root_2deriv_s(radius_a, arc_length, radius_b / 2)
430 << "][" << elliptic_root_2deriv_s(radius_a, arc_length, radius_b - radius_b * 0.5)
431 << "][" << elliptic_root_2deriv_s(radius_a, arc_length, radius_b - radius_b * 0.2)
432 << "][" << elliptic_root_2deriv_s(radius_a, arc_length, radius_b - radius_b * 0.1)
433 << "][" << elliptic_root_2deriv_s(radius_a, arc_length, radius_b - radius_b * 0.05)
434 << "][" << elliptic_root_2deriv_s(radius_a, arc_length, radius_b + radius_b * 0.05)
435 << "][" << elliptic_root_2deriv_s(radius_a, arc_length, radius_b + radius_b * 0.1)
436 << "][" << elliptic_root_2deriv_s(radius_a, arc_length, radius_b + radius_b * 0.2)
437 << "][" << elliptic_root_2deriv_s(radius_a, arc_length, radius_b + radius_b * 0.5)
438 << "][" << elliptic_root_2deriv_s(radius_a, arc_length, radius_b + radius_b)
439 << "][" << elliptic_root_2deriv_s(radius_a, arc_length, radius_b + radius_b * 5) << "]]\n]\n\n";
440
441 return boost::exit_success;
442 }
443 catch(std::exception ex)
444 {
445 std::cout << "exception thrown: " << ex.what() << std::endl;
446 return boost::exit_failure;
447 }
448 } // int main()
449
450