Comparison of Elliptic Integral Root Finding Algorithms

A second example compares four root finding algorithms for locating the second radius of an ellipse with first radius 28 and arc length 300, for four floating-point types, float, double, long double and a Boost.Multiprecision type cpp_bin_float_50.

Which is to say we're solving:

4xE(sqrt(1 - 28² / x²)) - 300 = 0

In each case the target accuracy was set using our "recommended" accuracy limits (or at least limits that make a good starting point - which is likely to give close to full accuracy without resorting to unnecessary iterations).

Function	Precision Requested
TOMS748	numeric_limits<T>::digits - 2
Newton	floor(numeric_limits<T>::digits * 0.6)
Halley	floor(numeric_limits<T>::digits * 0.4)
Schröder	floor(numeric_limits<T>::digits * 0.4)

Tests used Microsoft Visual Studio 2013 (Update 1) and GCC 4.9.1 using source code root_elliptic_finding.cpp.

The timing uncertainty (especially using MSVC) is at least 5% of normalized time 'Norm'.

To pick out the 'best' and 'worst' algorithms are highlighted in blue and red. More than one result can be 'best' when normalized times are indistinguishable within the uncertainty.

Program root_elliptic_finding.cpp, Microsoft Visual C++ version 14.1, Dinkumware standard library version 650, Win32 Compiled in optimise mode., _X86_SSE2

Table 10.12. root with radius 28 and arc length 300) for float, double, long double and cpp_bin_float_50 types, using _X86_SSE2

	float				double				long d				cpp50
Algo	Its	Times	Norm	Dis	Its	Times	Norm	Dis	Its	Times	Norm	Dis	Its	Times	Norm	Dis
TOMS748	6	906	2.07	0	9	1312	1.79	1	9	1281	1.75	1	11	1690625	1.52	-3
Newton	3	640	1.46	-1	4	875	1.19	1	4	843	1.15	1	5	1368750	1.23	0
Halley	2	437	1.00	0	3	734	1.00	3	3	734	1.00	3	4	1109375	1.00	0
Schröder	3	671	1.54	-1	6	1296	1.77	1	6	1406	1.92	1	5	1462500	1.32	-2

Program root_elliptic_finding.cpp, Microsoft Visual C++ version 12.0, Dinkumware standard library version 610, Win32 Compiled in optimise mode., _X64_AVX

Table 10.13. root with radius 28 and arc length 300) for float, double, long double and cpp_bin_float_50 types, using _X64_AVX

	float				double				long d				cpp50
Algo	Its	Times	Norm	Dis	Its	Times	Norm	Dis	Its	Times	Norm	Dis	Its	Times	Norm	Dis
TOMS748	5	500	1.33	-1	9	1046	1.72	1	9	1062	1.70	1	11	698437	1.54	-3
Newton	3	484	1.29	-1	4	734	1.21	1	4	687	1.10	1	5	545312	1.20	0
Halley	2	375	1.00	0	3	609	1.00	3	3	625	1.00	3	4	453125	1.00	0
Schröder	3	546	1.46	-1	6	1109	1.82	1	6	1187	1.90	1	5	564062	1.24	-2

Program root_elliptic_finding.cpp, GNU C++ version 7.1.0, GNU libstdc++ version 20170502, Win32 Compiled in optimise mode., _X64_SSE2

Table 10.14. root with radius 28 and arc length 300) for float, double, long double and cpp_bin_float_50 types, using _X64_SSE2

	float				double				long d				cpp50
Algo	Its	Times	Norm	Dis	Its	Times	Norm	Dis	Its	Times	Norm	Dis	Its	Times	Norm	Dis
TOMS748	5	328	1.24	-1	8	890	1.50	0	8	1234	1.61	4	11	487500	1.57	-3
Newton	3	359	1.35	-1	4	718	1.21	1	4	843	1.10	1	5	379687	1.22	0
Halley	2	265	1.00	0	3	593	1.00	1	3	765	1.00	7	4	310937	1.00	0
Schröder	3	343	1.29	-1	4	812	1.37	0	4	1046	1.37	3	5	390625	1.26	-2

Remarks:

The function being solved is now moderately expensive to call, and twice as expensive to call when obtaining the derivative than when not. Consequently there is very little improvement in moving from a derivative free method, to Newton iteration. However, once you've calculated the first derivative the second comes almost for free, consequently the third order methods (Halley) does much the best.
Of the two second order methods, Halley does best as would be expected: the Schroder method offers better guarantees of quadratic convergence, while Halley relies on a smooth function with a single root to give cubic convergence. It's not entirely clear why Schroder iteration often does worse than Newton.