]>
git.proxmox.com Git - ceph.git/blob - ceph/src/boost/libs/math/test/test_root_iterations.cpp
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)
8 #define BOOST_TEST_MAIN
9 #include <boost/test/unit_test.hpp>
10 #include <boost/test/tools/floating_point_comparison.hpp>
11 #include <boost/math/tools/roots.hpp>
12 #include <boost/test/results_collector.hpp>
13 #include <boost/test/unit_test.hpp>
14 #include <boost/math/special_functions/cbrt.hpp>
15 #include <boost/math/special_functions/beta.hpp>
19 #include "table_type.hpp"
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.
27 double operator ()( double x
)
29 double fx
= x
* x
* x
- a
; // Difference (estimate x^3 - a).
33 double a
; // to be 'cube_rooted'.
34 }; // template <class T> struct cbrt_functor_noderiv
36 // Using 1st derivative only Newton-Raphson
37 struct cbrt_functor_deriv
38 { // Functor also returning 1st derivative.
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.
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.
50 double a
; // to be 'cube_rooted'.
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,
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.
67 double a
; // to be 'cube_rooted'.
70 template < class T
, class Policy
>
71 struct ibeta_roots_1
// for first order algorithms
73 ibeta_roots_1 ( T _a
, T _b
, T t
, bool inv
= false )
74 : a ( _a
), b ( _b
), target ( t
), invert ( inv
) {}
76 T
operator ()( const T
& x
)
78 return boost :: math :: detail :: ibeta_imp ( a
, b
, x
, Policy (), invert
, true ) - target
;
85 template < class T
, class Policy
>
86 struct ibeta_roots_2
// for second order algorithms
88 ibeta_roots_2 ( T _a
, T _b
, T t
, bool inv
= false )
89 : a ( _a
), b ( _b
), target ( t
), invert ( inv
) {}
91 boost :: math :: tuple
< T
, T
> operator ()( const T
& x
)
93 typedef boost :: math :: lanczos :: lanczos
< T
, Policy
> S
;
94 typedef typename
S :: type L
;
95 T f
= boost :: math :: detail :: ibeta_imp ( a
, b
, x
, Policy (), invert
, true ) - target
;
97 - boost :: math :: detail :: ibeta_power_terms ( b
, a
, 1 - x
, x
, L (), true , Policy ())
98 : boost :: math :: detail :: ibeta_power_terms ( a
, b
, x
, 1 - x
, L (), true , Policy ());
101 y
= boost :: math :: tools :: min_value
< T
>() * 8 ;
104 // make sure we don't have a zero derivative:
106 f1
= ( invert
? - 1 : 1 ) * boost :: math :: tools :: min_value
< T
>() * 64 ;
108 return boost :: math :: make_tuple ( f
, f1
);
115 template < class T
, class Policy
>
116 struct ibeta_roots_3
// for third order algorithms
118 ibeta_roots_3 ( T _a
, T _b
, T t
, bool inv
= false )
119 : a ( _a
), b ( _b
), target ( t
), invert ( inv
) {}
121 boost :: math :: tuple
< T
, T
, T
> operator ()( const T
& x
)
123 typedef typename
boost :: math :: lanczos :: lanczos
< T
, Policy
>:: type L
;
124 T f
= boost :: math :: detail :: ibeta_imp ( a
, b
, x
, Policy (), invert
, true ) - target
;
126 - boost :: math :: detail :: ibeta_power_terms ( b
, a
, 1 - x
, x
, L (), true , Policy ())
127 : boost :: math :: detail :: ibeta_power_terms ( a
, b
, x
, 1 - x
, L (), true , Policy ());
130 y
= boost :: math :: tools :: min_value
< T
>() * 8 ;
132 T f2
= f1
* (- y
* a
+ ( b
- 2 ) * x
+ 1 ) / ( y
* x
);
136 // make sure we don't have a zero derivative:
138 f1
= ( invert
? - 1 : 1 ) * boost :: math :: tools :: min_value
< T
>() * 64 ;
140 return boost :: math :: make_tuple ( f
, f1
, f2
);
148 BOOST_AUTO_TEST_CASE ( test_main
)
150 int newton_limits
= static_cast < int >( std :: numeric_limits
< double >:: digits
* 0.6 );
153 std :: uintmax_t iters
;
159 double result
= boost :: math :: cbrt ( arg
);
161 // Start with a really bad guess 5 times below the result:
166 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
);
167 BOOST_CHECK_CLOSE_FRACTION (( r
. first
+ r
. second
) / 2 , result
, std :: numeric_limits
< double >:: epsilon () * 4 );
168 BOOST_CHECK_LE ( iters
, 14 );
171 dr
= boost :: math :: tools :: newton_raphson_iterate ( cbrt_functor_deriv ( arg
), guess
, guess
/ 2 , result
* 10 , newton_limits
, iters
);
172 BOOST_CHECK_CLOSE_FRACTION ( dr
, result
, std :: numeric_limits
< double >:: epsilon () * 2 );
173 BOOST_CHECK_LE ( iters
, 12 );
176 dr
= boost :: math :: tools :: halley_iterate ( cbrt_functor_2deriv ( arg
), guess
, result
/ 10 , result
* 10 , newton_limits
, iters
);
177 BOOST_CHECK_CLOSE_FRACTION ( dr
, result
, std :: numeric_limits
< double >:: epsilon () * 2 );
178 BOOST_CHECK_LE ( iters
, 7 );
181 dr
= boost :: math :: tools :: schroder_iterate ( cbrt_functor_2deriv ( arg
), guess
, result
/ 10 , result
* 10 , newton_limits
, iters
);
182 BOOST_CHECK_CLOSE_FRACTION ( dr
, result
, std :: numeric_limits
< double >:: epsilon () * 2 );
183 BOOST_CHECK_LE ( iters
, 11 );
185 // Over again with a bad guess 5 times larger than the result:
189 r
= boost :: math :: tools :: bracket_and_solve_root ( cbrt_functor_noderiv ( arg
), guess
, 2.0 , true , boost :: math :: tools :: eps_tolerance
< double >(), iters
);
190 BOOST_CHECK_CLOSE_FRACTION (( r
. first
+ r
. second
) / 2 , result
, std :: numeric_limits
< double >:: epsilon () * 4 );
191 BOOST_CHECK_LE ( iters
, 14 );
194 dr
= boost :: math :: tools :: newton_raphson_iterate ( cbrt_functor_deriv ( arg
), guess
, result
/ 10 , result
* 10 , newton_limits
, iters
);
195 BOOST_CHECK_CLOSE_FRACTION ( dr
, result
, std :: numeric_limits
< double >:: epsilon () * 2 );
196 BOOST_CHECK_LE ( iters
, 12 );
199 dr
= boost :: math :: tools :: halley_iterate ( cbrt_functor_2deriv ( arg
), guess
, result
/ 10 , result
* 10 , newton_limits
, iters
);
200 BOOST_CHECK_CLOSE_FRACTION ( dr
, result
, std :: numeric_limits
< double >:: epsilon () * 2 );
201 BOOST_CHECK_LE ( iters
, 7 );
204 dr
= boost :: math :: tools :: schroder_iterate ( cbrt_functor_2deriv ( arg
), guess
, result
/ 10 , result
* 10 , newton_limits
, iters
);
205 BOOST_CHECK_CLOSE_FRACTION ( dr
, result
, std :: numeric_limits
< double >:: epsilon () * 2 );
206 BOOST_CHECK_LE ( iters
, 11 );
208 // A much better guess, 1% below result:
211 guess
= result
* 0.9 ;
212 r
= boost :: math :: tools :: bracket_and_solve_root ( cbrt_functor_noderiv ( arg
), guess
, 2.0 , true , boost :: math :: tools :: eps_tolerance
< double >(), iters
);
213 BOOST_CHECK_CLOSE_FRACTION (( r
. first
+ r
. second
) / 2 , result
, std :: numeric_limits
< double >:: epsilon () * 4 );
214 BOOST_CHECK_LE ( iters
, 12 );
217 dr
= boost :: math :: tools :: newton_raphson_iterate ( cbrt_functor_deriv ( arg
), guess
, result
/ 10 , result
* 10 , newton_limits
, iters
);
218 BOOST_CHECK_CLOSE_FRACTION ( dr
, result
, std :: numeric_limits
< double >:: epsilon () * 2 );
219 BOOST_CHECK_LE ( iters
, 5 );
222 dr
= boost :: math :: tools :: halley_iterate ( cbrt_functor_2deriv ( arg
), guess
, result
/ 10 , result
* 10 , newton_limits
, iters
);
223 BOOST_CHECK_CLOSE_FRACTION ( dr
, result
, std :: numeric_limits
< double >:: epsilon () * 2 );
224 BOOST_CHECK_LE ( iters
, 3 );
227 dr
= boost :: math :: tools :: schroder_iterate ( cbrt_functor_2deriv ( arg
), guess
, result
/ 10 , result
* 10 , newton_limits
, iters
);
228 BOOST_CHECK_CLOSE_FRACTION ( dr
, result
, std :: numeric_limits
< double >:: epsilon () * 2 );
229 BOOST_CHECK_LE ( iters
, 4 );
231 // A much better guess, 1% above result:
234 guess
= result
* 1.1 ;
235 r
= boost :: math :: tools :: bracket_and_solve_root ( cbrt_functor_noderiv ( arg
), guess
, 2.0 , true , boost :: math :: tools :: eps_tolerance
< double >(), iters
);
236 BOOST_CHECK_CLOSE_FRACTION (( r
. first
+ r
. second
) / 2 , result
, std :: numeric_limits
< double >:: epsilon () * 4 );
237 BOOST_CHECK_LE ( iters
, 12 );
240 dr
= boost :: math :: tools :: newton_raphson_iterate ( cbrt_functor_deriv ( arg
), guess
, result
/ 10 , result
* 10 , newton_limits
, iters
);
241 BOOST_CHECK_CLOSE_FRACTION ( dr
, result
, std :: numeric_limits
< double >:: epsilon () * 2 );
242 BOOST_CHECK_LE ( iters
, 5 );
245 dr
= boost :: math :: tools :: halley_iterate ( cbrt_functor_2deriv ( arg
), guess
, result
/ 10 , result
* 10 , newton_limits
, iters
);
246 BOOST_CHECK_CLOSE_FRACTION ( dr
, result
, std :: numeric_limits
< double >:: epsilon () * 2 );
247 BOOST_CHECK_LE ( iters
, 3 );
250 dr
= boost :: math :: tools :: schroder_iterate ( cbrt_functor_2deriv ( arg
), guess
, result
/ 10 , result
* 10 , newton_limits
, iters
);
251 BOOST_CHECK_CLOSE_FRACTION ( dr
, result
, std :: numeric_limits
< double >:: epsilon () * 2 );
252 BOOST_CHECK_LE ( iters
, 4 );
258 // Test ibeta as this triggers all the pathological cases!
265 # include "ibeta_small_data.ipp"
267 for ( unsigned i
= 0 ; i
< ibeta_small_data
. size (); ++ i
)
270 // These inverse tests are thrown off if the output of the
271 // incomplete beta is too close to 1: basically there is insuffient
272 // information left in the value we're using as input to the inverse
273 // to be able to get back to the original value.
275 if ( ibeta_small_data
[ i
][ 5 ] == 0 )
278 dr
= boost :: math :: tools :: newton_raphson_iterate ( ibeta_roots_2
< double , boost :: math :: policies :: policy
<> >( ibeta_small_data
[ i
][ 0 ], ibeta_small_data
[ i
][ 1 ], ibeta_small_data
[ i
][ 5 ]), 0.5 , 0.0 , 1.0 , 53 , iters
);
279 BOOST_CHECK_EQUAL ( dr
, 0.0 );
280 BOOST_CHECK_LE ( iters
, 27 );
282 dr
= boost :: math :: tools :: halley_iterate ( ibeta_roots_3
< double , boost :: math :: policies :: policy
<> >( ibeta_small_data
[ i
][ 0 ], ibeta_small_data
[ i
][ 1 ], ibeta_small_data
[ i
][ 5 ]), 0.5 , 0.0 , 1.0 , 53 , iters
);
283 BOOST_CHECK_EQUAL ( dr
, 0.0 );
284 BOOST_CHECK_LE ( iters
, 10 );
286 else if (( 1 - ibeta_small_data
[ i
][ 5 ] > 0.001 )
287 && ( fabs ( ibeta_small_data
[ i
][ 5 ]) > 2 * boost :: math :: tools :: min_value
< double >()))
290 double result
= ibeta_small_data
[ i
][ 2 ];
291 dr
= boost :: math :: tools :: newton_raphson_iterate ( ibeta_roots_2
< double , boost :: math :: policies :: policy
<> >( ibeta_small_data
[ i
][ 0 ], ibeta_small_data
[ i
][ 1 ], ibeta_small_data
[ i
][ 5 ]), 0.5 , 0.0 , 1.0 , 53 , iters
);
292 BOOST_CHECK_CLOSE_FRACTION ( dr
, result
, std :: numeric_limits
< double >:: epsilon () * 200 );
293 #if defined(BOOST_MSVC) && (BOOST_MSVC == 1600)
294 BOOST_CHECK_LE ( iters
, 40 );
296 BOOST_CHECK_LE ( iters
, 27 );
299 result
= ibeta_small_data
[ i
][ 2 ];
300 dr
= boost :: math :: tools :: halley_iterate ( ibeta_roots_3
< double , boost :: math :: policies :: policy
<> >( ibeta_small_data
[ i
][ 0 ], ibeta_small_data
[ i
][ 1 ], ibeta_small_data
[ i
][ 5 ]), 0.5 , 0.0 , 1.0 , 53 , iters
);
301 BOOST_CHECK_CLOSE_FRACTION ( dr
, result
, std :: numeric_limits
< double >:: epsilon () * 200 );
302 #if defined(__PPC__) || defined(__aarch64__)
303 BOOST_CHECK_LE ( iters
, 55 );
305 BOOST_CHECK_LE ( iters
, 40 );
308 else if ( 1 == ibeta_small_data
[ i
][ 5 ])
311 dr
= boost :: math :: tools :: newton_raphson_iterate ( ibeta_roots_2
< double , boost :: math :: policies :: policy
<> >( ibeta_small_data
[ i
][ 0 ], ibeta_small_data
[ i
][ 1 ], ibeta_small_data
[ i
][ 5 ]), 0.5 , 0.0 , 1.0 , 53 , iters
);
312 BOOST_CHECK_EQUAL ( dr
, 1.0 );
313 BOOST_CHECK_LE ( iters
, 27 );
315 dr
= boost :: math :: tools :: halley_iterate ( ibeta_roots_3
< double , boost :: math :: policies :: policy
<> >( ibeta_small_data
[ i
][ 0 ], ibeta_small_data
[ i
][ 1 ], ibeta_small_data
[ i
][ 5 ]), 0.5 , 0.0 , 1.0 , 53 , iters
);
316 BOOST_CHECK_EQUAL ( dr
, 1.0 );
317 BOOST_CHECK_LE ( iters
, 10 );