1[/ 2Copyright 2015 Paul A. Bristow. 3Copyright 2015 John Maddock. 4Distributed under the Boost Software License, Version 1.0. 5(See accompanying file LICENSE_1_0.txt or copy at 6http://www.boost.org/LICENSE_1_0.txt). 7] 8 9[section:root_comparison Comparison of Root Finding Algorithms] 10 11[section:cbrt_comparison Comparison of Cube Root Finding Algorithms] 12 13In the table below, the cube root of 28 was computed for three __fundamental_types floating-point types, 14and one __multiprecision type __cpp_bin_float using 50 decimal digit precision, using four algorithms. 15 16The 'exact' answer was computed using a 100 decimal digit type: 17 18 cpp_bin_float_100 full_answer ("3.036588971875662519420809578505669635581453977248111123242141654169177268411884961770250390838097895"); 19 20Times were measured using __boost_timer using `class cpu_timer`. 21 22* ['Its] is the number of iterations taken to find the root. 23* ['Times] is the CPU time-taken in arbitrary units. 24* ['Norm] is a normalized time, in comparison to the quickest algorithm (with value 1.00). 25* ['Dis] is the distance from the nearest representation of the 'exact' root in bits. 26Distance from the 'exact' answer is measured by using function __float_distance. 27One or two bits distance means that all results are effectively 'correct'. 28Zero means 'exact' - the nearest __representable value for the floating-point type. 29 30The cube-root function is a simple function, and is a contrived example for root-finding. 31It does allow us to investigate some of the factors controlling efficiency that may be extrapolated to 32more complex functions. 33 34The program used was [@ ../../example/root_finding_algorithms.cpp root_finding_algorithms.cpp]. 35100000 evaluations of each floating-point type and algorithm were used and the CPU times were 36judged from repeat runs to have an uncertainty of 10 %. Comparing MSVC for `double` and `long double` 37(which are identical on this platform) may give a guide to uncertainty of timing. 38 39The requested precision was set as follows: 40 41[table 42[[Function][Precision Requested]] 43[[TOMS748][numeric_limits<T>::digits - 2]] 44[[Newton][floor(numeric_limits<T>::digits * 0.6)]] 45[[Halley][floor(numeric_limits<T>::digits * 0.4)]] 46[[Schr'''ö'''der][floor(numeric_limits<T>::digits * 0.4)]] 47] 48 49* The C++ Standard cube root function [@http://en.cppreference.com/w/cpp/numeric/math/cbrt std::cbrt] 50is only defined for built-in or fundamental types, 51so cannot be used with any User-Defined floating-point types like __multiprecision. 52This, and that the cube function is so impeccably-behaved, 53allows the implementer to use many tricks to achieve a fast computation. 54On some platforms,`std::cbrt` appeared several times as quick as the more general `boost::math::cbrt`, 55on other platforms / compiler options `boost::math::cbrt` is noticeably faster. In general, the results are highly 56dependent on the code-generation / processor architecture selection compiler options used. One can 57assume that the standard library will have been compiled with options ['nearly] optimal for the platform 58it was installed on, where as the user has more choice over the options used for Boost.Math. Pick something 59too general/conservative and performance suffers, while selecting options that make use of the latest 60instruction set opcodes speed's things up noticeably. 61 62* Two compilers in optimise mode were compared: GCC 4.9.1 using Netbeans IDS 63and Microsoft Visual Studio 2013 (Update 1) on the same hardware. 64The number of iterations seemed consistent, but the relative run-times surprisingly different. 65 66* `boost::math::cbrt` allows use with ['any user-defined floating-point type], conveniently 67__multiprecision. It too can take some advantage of the good-behaviour of the cube function, 68compared to the more general implementation in the nth root-finding examples. For example, 69it uses a polynomial approximation to generate a better guess than dividing the exponent by three, 70and can avoid the complex checks in __newton required to prevent the 71search going wildly off-track. For a known precision, it may also be possible to 72fix the number of iterations, allowing inlining and loop unrolling. It also 73algebraically simplifies the Halley steps leading to a big reduction in the 74number of floating point operations required compared to a "black box" implementation 75that calculates the derivatives separately and then combines them in the Halley code. 76Typically, it was found that computation using type `double` 77took a few times longer when using the various root-finding algorithms directly rather 78than the hand coded/optimized `cbrt` routine. 79 80* The importance of getting a good guess can be seen by the iteration count for the multiprecision case: 81here we "cheat" a little and use the cube-root calculated to double precision as the initial guess. 82The limitation of this tactic is that the range of possible (exponent) values may be less than the multiprecision type. 83 84* For __fundamental_types, there was little to choose between the three derivative methods, 85but for __cpp_bin_float, __newton was twice as fast. Note that the cube-root is an extreme 86test case as the cost of calling the functor is so cheap that the runtimes are largely 87dominated by the complexity of the iteration code. 88 89* Compiling with optimisation halved computation times, and any differences between algorithms 90became nearly negligible. The optimisation speed-up of the __TOMS748 was especially noticeable. 91 92* Using a multiprecision type like `cpp_bin_float_50` for a precision of 50 decimal digits 93took a lot longer, as expected because most computation 94uses software rather than 64-bit floating-point hardware. 95Speeds are often more than 50 times slower. 96 97* Using `cpp_bin_float_50`, __TOMS748 was much slower showing the benefit of using derivatives. 98__newton was found to be twice as quick as either of the second-derivative methods: 99this is an extreme case though, the function and its derivatives are so cheap to compute that we're 100really measuring the complexity of the boilerplate root-finding code. 101 102* For multiprecision types only one or two extra ['iterations] are needed to get the remaining 35 digits, whatever the algorithm used. 103(The time taken was of course much greater for these types). 104 105* Using a 100 decimal-digit type only doubled the time and required only a very few more iterations, 106so the cost of extra precision is mainly the underlying cost of computing more digits, 107not in the way the algorithm works. This confirms previous observations using __NTL high-precision types. 108 109[include root_comparison_tables_msvc.qbk] 110[include root_comparison_tables_gcc.qbk] 111 112[endsect] [/section:cbrt_comparison Comparison of Cube Root Finding Algorithms] 113 114[section:root_n_comparison Comparison of Nth-root Finding Algorithms] 115 116A second example compares four generalized nth-root finding algorithms for various n-th roots (5, 7 and 13) 117of a single value 28.0, for four floating-point types, `float`, `double`, 118`long double` and a __multiprecision type `cpp_bin_float_50`. 119In each case the target accuracy was set using our "recommended" accuracy limits 120(or at least limits that make a good starting point - which is likely to give 121close to full accuracy without resorting to unnecessary iterations). 122 123[table 124[[Function][Precision Requested]] 125[[TOMS748][numeric_limits<T>::digits - 2]] 126[[Newton][floor(numeric_limits<T>::digits * 0.6)]] 127[[Halley][floor(numeric_limits<T>::digits * 0.4)]] 128[[Schr'''ö'''der][floor(numeric_limits<T>::digits * 0.4)]] 129] 130Tests used Microsoft Visual Studio 2013 (Update 1) and GCC 4.9.1 using source code 131[@../../example/root_n_finding_algorithms.cpp root_n_finding_algorithms.cpp]. 132 133The timing uncertainty (especially using MSVC) is at least 5% of normalized time 'Norm'. 134 135To pick out the 'best' and 'worst' algorithms are highlighted in blue and red. 136More than one result can be 'best' when normalized times are indistinguishable 137within the uncertainty. 138 139[/include roots_table_100_msvc.qbk] 140[/include roots_table_75_msvc.qbk] 141 142[/include roots_table_75_msvc_X86.qbk] 143[/include roots_table_100_msvc_X86.qbk] 144 145[/include roots_table_100_msvc_AVX.qbk] 146[/include roots_table_75_msvc_AVX.qbk] 147 148[/include roots_table_75_msvc_X86_SSE2.qbk] 149[/include roots_table_100_msvc_X86_SSE2.qbk] 150 151[/include roots_table_100_gcc_X64_SSE2.qbk] 152[/include roots_table_75_gcc_X64_SSE2.qbk] 153 154[/include type_info_table_100_msvc.qbk] 155[/include type_info_table_75_msvc.qbk] 156 157[include roots_table_100_msvc_X86_SSE2.qbk] 158[include roots_table_100_msvc_X64_AVX.qbk] 159[include roots_table_100_gcc_X64_SSE2.qbk] 160 161Some tentative conclusions can be drawn from this limited exercise. 162 163* Perhaps surprisingly, there is little difference between the various algorithms for __fundamental_types floating-point types. 164Using the first derivatives (__newton) is usually the best, but while the improvement over the no-derivative 165__TOMS748 is considerable in number of iterations, but little in execution time. This reflects the fact that the function 166we are finding the root for is trivial to evaluate, so runtimetimes are dominated by the time taken by the boilerplate code 167in each method. 168 169* The extra cost of evaluating the second derivatives (__halley or __schroder) is usually too much for any net benefit: 170as with the cube root, these functors are so cheap to evaluate that the runtime is largely dominated by the 171complexity of the root finding method. 172 173* For a __multiprecision floating-point type, the __newton is a clear winner with a several-fold gain over __TOMS748, 174and again no improvement from the second-derivative algorithms. 175 176* The run-time of 50 decimal-digit __multiprecision is about 30-fold greater than `double`. 177 178* The column 'dis' showing the number of bits distance from the correct result. 179The Newton-Raphson algorithm shows a bit or two better accuracy than __TOMS748. 180 181* The goodness of the 'guess' is especially crucial for __multiprecision. 182Separate experiments show that evaluating the 'guess' using `double` allows 183convergence to the final exact result in one or two iterations. 184So in this contrived example, crudely dividing the exponent by N for a 'guess', 185it would be far better to use a `pow<double>` or , 186if more precise `pow<long double>`, function to estimate a 'guess'. 187The limitation of this tactic is that the range of possible (exponent) values may be less than the multiprecision type. 188 189* Using floating-point extension __SSE2 made a modest ten-percent speedup. 190 191*Using MSVC, there was some improvement using 64-bit, markedly for __multiprecision. 192 193* The GCC compiler 4.9.1 using 64-bit was at least five-folder faster that 32-bit, 194apparently reflecting better optimization. 195 196Clearly, your mileage [*will vary], but in summary, __newton seems the first choice of algorithm, 197and effort to find a good 'guess' the first speed-up target, especially for __multiprecision. 198And of course, compiler optimisation is crucial for speed. 199 200[endsect] [/section:root_n_comparison Comparison of Nth-root Finding Algorithms] 201 202[section:elliptic_comparison Comparison of Elliptic Integral Root Finding Algorithms] 203 204A second example compares four root finding algorithms for locating 205the second radius of an ellipse with first radius 28 and arc length 300, 206for four floating-point types, `float`, `double`, 207`long double` and a __multiprecision type `cpp_bin_float_50`. 208 209Which is to say we're solving: 210 211[pre 4xE(sqrt(1 - 28[super 2] / x[super 2])) - 300 = 0] 212 213In each case the target accuracy was set using our "recommended" accuracy limits 214(or at least limits that make a good starting point - which is likely to give 215close to full accuracy without resorting to unnecessary iterations). 216 217[table 218[[Function][Precision Requested]] 219[[TOMS748][numeric_limits<T>::digits - 2]] 220[[Newton][floor(numeric_limits<T>::digits * 0.6)]] 221[[Halley][floor(numeric_limits<T>::digits * 0.4)]] 222[[Schr'''ö'''der][floor(numeric_limits<T>::digits * 0.4)]] 223] 224Tests used Microsoft Visual Studio 2013 (Update 1) and GCC 4.9.1 using source code 225[@../../example/root_elliptic_finding.cpp root_elliptic_finding.cpp]. 226 227The timing uncertainty (especially using MSVC) is at least 5% of normalized time 'Norm'. 228 229To pick out the 'best' and 'worst' algorithms are highlighted in blue and red. 230More than one result can be 'best' when normalized times are indistinguishable 231within the uncertainty. 232 233[include elliptic_table_100_msvc_X86_SSE2.qbk] 234[include elliptic_table_100_msvc_X64_AVX.qbk] 235[include elliptic_table_100_gcc_X64_SSE2.qbk] 236 237Remarks: 238 239* The function being solved is now moderately expensive to call, and twice as expensive to call 240when obtaining the derivative than when not. Consequently there is very little improvement in moving 241from a derivative free method, to Newton iteration. However, once you've calculated the first derivative 242the second comes almost for free, consequently the third order methods (Halley) does much the best. 243* Of the two second order methods, Halley does best as would be expected: the Schroder method offers better 244guarantees of ['quadratic] convergence, while Halley relies on a smooth function with a single root to 245give ['cubic] convergence. It's not entirely clear why Schroder iteration often does worse than Newton. 246 247[endsect] [/section:elliptic_comparison Comparison of Elliptic Integral Root Finding Algorithms] 248 249[endsect] [/section:root_comparison Comparison of Root Finding Algorithms] 250