1 // (C) Copyright John Maddock 2006.
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)
6 #define BOOST_TEST_MODULE foobar
7 #define BOOST_UBLAS_TYPE_CHECK_EPSILON (type_traits<real_type>::type_sqrt (boost::math::tools::epsilon <real_type>()))
8 #define BOOST_UBLAS_TYPE_CHECK_MIN (type_traits<real_type>::type_sqrt ( boost::math::tools::min_value<real_type>()))
9 #define BOOST_UBLAS_NDEBUG
11 #include "multiprecision.hpp"
13 #include <boost/math/tools/remez.hpp>
14 #include <boost/math/tools/test.hpp>
15 #include <boost/math/special_functions/binomial.hpp>
16 #include <boost/spirit/include/classic_core.hpp>
17 #include <boost/spirit/include/classic_actor.hpp>
18 #include <boost/lexical_cast.hpp>
22 #include <boost/test/included/unit_test.hpp> // for test_main
23 #include <boost/multiprecision/cpp_bin_float.hpp>
26 extern mp_type
f(const mp_type
& x
, int variant
);
27 extern void show_extra(
28 const boost::math::tools::polynomial
<mp_type
>& n
,
29 const boost::math::tools::polynomial
<mp_type
>& d
,
30 const mp_type
& x_offset
,
31 const mp_type
& y_offset
,
34 using namespace boost::spirit::classic
;
36 mp_type
a(0), b(1); // range to optimise over
41 int target_precision
= boost::math::tools::digits
<long double>();
42 int working_precision
= target_precision
* 2;
47 mp_type
x_offset(0), y_offset(0), x_scale(1);
50 boost::shared_ptr
<boost::math::tools::remez_minimax
<mp_type
> > p_remez
;
52 mp_type
the_function(const mp_type
& val
)
54 return f(x_scale
* (val
+ x_offset
), variant
) + y_offset
;
57 void step_some(unsigned count
)
60 set_working_precision(working_precision
);
64 // If we have an automatic y-offset calculate it now:
69 fa
= f(x_scale
* (a
+ x_offset
), variant
);
70 fb
= f(x_scale
* (b
+ x_offset
), variant
);
71 fm
= f(x_scale
* ((a
+b
)/2 + x_offset
), variant
);
72 y_offset
= -(fa
+ fb
+ fm
) / 3;
73 set_output_precision(5);
74 std::cout
<< "Setting auto-y-offset to " << y_offset
<< std::endl
;
77 // Truncate offsets to float precision:
79 x_offset
= round_to_precision(x_offset
, 20);
80 y_offset
= round_to_precision(y_offset
, 20);
82 // Construct new Remez state machine:
84 p_remez
.reset(new boost::math::tools::remez_minimax
<mp_type
>(
92 std::cout
<< "Max error in interpolated form: " << std::setprecision(3) << std::scientific
<< boost::math::tools::real_cast
<double>(p_remez
->max_error()) << std::endl
;
94 // Signal that we've started:
99 for(i
= 0; i
< count
; ++i
)
101 std::cout
<< "Stepping..." << std::endl
;
102 p_remez
->set_brake(brake
);
103 mp_type r
= p_remez
->iterate();
104 set_output_precision(3);
106 << "Maximum Deviation Found: " << std::setprecision(3) << std::scientific
<< boost::math::tools::real_cast
<double>(p_remez
->max_error()) << std::endl
107 << "Expected Error Term: " << std::setprecision(3) << std::scientific
<< boost::math::tools::real_cast
<double>(p_remez
->error_term()) << std::endl
108 << "Maximum Relative Change in Control Points: " << std::setprecision(3) << std::scientific
<< boost::math::tools::real_cast
<double>(r
) << std::endl
;
111 catch(const std::exception
& e
)
113 std::cout
<< "Step failed with exception: " << e
.what() << std::endl
;
117 void step(const char*, const char*)
122 void show(const char*, const char*)
124 set_working_precision(working_precision
);
127 boost::math::tools::polynomial
<mp_type
> n
= p_remez
->numerator();
128 boost::math::tools::polynomial
<mp_type
> d
= p_remez
->denominator();
129 std::vector
<mp_type
> cn
= n
.chebyshev();
130 std::vector
<mp_type
> cd
= d
.chebyshev();
131 int prec
= 2 + (target_precision
* 3010LL)/10000;
132 std::cout
<< std::scientific
<< std::setprecision(prec
);
133 set_output_precision(prec
);
134 boost::numeric::ublas::vector
<mp_type
> v
= p_remez
->zero_points();
136 std::cout
<< " Zeros = {\n";
138 for(i
= 0; i
< v
.size(); ++i
)
140 std::cout
<< " " << v
[i
] << std::endl
;
144 v
= p_remez
->chebyshev_points();
145 std::cout
<< " Chebeshev Control Points = {\n";
146 for(i
= 0; i
< v
.size(); ++i
)
148 std::cout
<< " " << v
[i
] << std::endl
;
152 std::cout
<< "X offset: " << x_offset
<< std::endl
;
153 std::cout
<< "X scale: " << x_scale
<< std::endl
;
154 std::cout
<< "Y offset: " << y_offset
<< std::endl
;
156 std::cout
<< "P = {";
157 for(i
= 0; i
< n
.size(); ++i
)
159 std::cout
<< " " << n
[i
] << "L," << std::endl
;
163 std::cout
<< "Q = {";
164 for(i
= 0; i
< d
.size(); ++i
)
166 std::cout
<< " " << d
[i
] << "L," << std::endl
;
170 std::cout
<< "CP = {";
171 for(i
= 0; i
< cn
.size(); ++i
)
173 std::cout
<< " " << cn
[i
] << "L," << std::endl
;
177 std::cout
<< "CQ = {";
178 for(i
= 0; i
< cd
.size(); ++i
)
180 std::cout
<< " " << cd
[i
] << "L," << std::endl
;
184 show_extra(n
, d
, x_offset
, y_offset
, variant
);
188 std::cerr
<< "Nothing to display" << std::endl
;
192 void do_graph(unsigned points
)
194 set_working_precision(working_precision
);
195 mp_type step
= (b
- a
) / (points
- 1);
199 set_output_precision(10);
200 std::cout
<< std::setprecision(10) << std::setw(30) << std::left
201 << boost::lexical_cast
<std::string
>(x
) << the_function(x
) << std::endl
;
205 std::cout
<< std::setprecision(10) << std::setw(30) << std::left
206 << boost::lexical_cast
<std::string
>(b
) << the_function(b
) << std::endl
;
209 void graph(const char*, const char*)
215 mp_type
convert_to_rr(const T
& val
)
219 template <class Backend
, boost::multiprecision::expression_template_option ET
>
220 mp_type
convert_to_rr(const boost::multiprecision::number
<Backend
, ET
>& val
)
222 return boost::lexical_cast
<mp_type
>(val
.str());
226 void do_test(T
, const char* name
)
228 set_working_precision(working_precision
);
232 // We want to test the approximation at fixed precision:
233 // either float, double or long double. Begin by getting the
236 boost::math::tools::polynomial
<T
> n
, d
;
237 boost::math::tools::polynomial
<mp_type
> nr
, dr
;
238 nr
= p_remez
->numerator();
239 dr
= p_remez
->denominator();
243 std::vector
<mp_type
> cn1
, cd1
;
244 cn1
= nr
.chebyshev();
245 cd1
= dr
.chebyshev();
246 std::vector
<T
> cn
, cd
;
247 for(unsigned i
= 0; i
< cn1
.size(); ++i
)
249 cn
.push_back(boost::math::tools::real_cast
<T
>(cn1
[i
]));
251 for(unsigned i
= 0; i
< cd1
.size(); ++i
)
253 cd
.push_back(boost::math::tools::real_cast
<T
>(cd1
[i
]));
256 // We'll test at the Chebeshev control points which is where
257 // (in theory) the largest deviation should occur. For good
258 // measure we'll test at the zeros as well:
260 boost::numeric::ublas::vector
<mp_type
>
261 zeros(p_remez
->zero_points()),
262 cheb(p_remez
->chebyshev_points());
264 mp_type
max_error(0), cheb_max_error(0);
267 // Do the tests at the zeros:
269 std::cout
<< "Starting tests at " << name
<< " precision...\n";
270 std::cout
<< "Absissa Error (Poly) Error (Cheb)\n";
271 for(unsigned i
= 0; i
< zeros
.size(); ++i
)
273 mp_type true_result
= the_function(zeros
[i
]);
274 T absissa
= boost::math::tools::real_cast
<T
>(zeros
[i
]);
275 mp_type test_result
= convert_to_rr(n
.evaluate(absissa
) / d
.evaluate(absissa
));
276 mp_type cheb_result
= convert_to_rr(boost::math::tools::evaluate_chebyshev(cn
, absissa
) / boost::math::tools::evaluate_chebyshev(cd
, absissa
));
277 mp_type err
, cheb_err
;
280 err
= boost::math::tools::relative_error(test_result
, true_result
);
281 cheb_err
= boost::math::tools::relative_error(cheb_result
, true_result
);
285 err
= fabs(test_result
- true_result
);
286 cheb_err
= fabs(cheb_result
- true_result
);
290 if(cheb_err
> cheb_max_error
)
291 cheb_max_error
= cheb_err
;
292 std::cout
<< std::setprecision(6) << std::setw(15) << std::left
<< absissa
293 << std::setw(15) << std::left
<< boost::math::tools::real_cast
<T
>(err
) << boost::math::tools::real_cast
<T
>(cheb_err
) << std::endl
;
296 // Do the tests at the Chebeshev control points:
298 for(unsigned i
= 0; i
< cheb
.size(); ++i
)
300 mp_type true_result
= the_function(cheb
[i
]);
301 T absissa
= boost::math::tools::real_cast
<T
>(cheb
[i
]);
302 mp_type test_result
= convert_to_rr(n
.evaluate(absissa
) / d
.evaluate(absissa
));
303 mp_type cheb_result
= convert_to_rr(boost::math::tools::evaluate_chebyshev(cn
, absissa
) / boost::math::tools::evaluate_chebyshev(cd
, absissa
));
304 mp_type err
, cheb_err
;
307 err
= boost::math::tools::relative_error(test_result
, true_result
);
308 cheb_err
= boost::math::tools::relative_error(cheb_result
, true_result
);
312 err
= fabs(test_result
- true_result
);
313 cheb_err
= fabs(cheb_result
- true_result
);
317 std::cout
<< std::setprecision(6) << std::setw(15) << std::left
<< absissa
318 << std::setw(15) << std::left
<< boost::math::tools::real_cast
<T
>(err
) <<
319 boost::math::tools::real_cast
<T
>(cheb_err
) << std::endl
;
321 std::string msg
= "Max Error found at ";
323 msg
+= " precision = ";
324 msg
.append(62 - 17 - msg
.size(), ' ');
325 std::cout
<< msg
<< std::setprecision(6) << "Poly: " << std::setw(20) << std::left
326 << boost::math::tools::real_cast
<T
>(max_error
) << "Cheb: " << boost::math::tools::real_cast
<T
>(cheb_max_error
) << std::endl
;
330 std::cout
<< "Nothing to test: try converging an approximation first!!!" << std::endl
;
334 void test_float(const char*, const char*)
336 do_test(float(0), "float");
339 void test_double(const char*, const char*)
341 do_test(double(0), "double");
344 void test_long(const char*, const char*)
346 do_test((long double)(0), "long double");
349 void test_float80(const char*, const char*)
351 do_test((boost::multiprecision::cpp_bin_float_double_extended
)(0), "float80");
354 void test_float128(const char*, const char*)
356 do_test((boost::multiprecision::cpp_bin_float_quad
)(0), "float128");
359 void test_all(const char*, const char*)
361 do_test(float(0), "float");
362 do_test(double(0), "double");
363 do_test((long double)(0), "long double");
367 void do_test_n(T
, const char* name
, unsigned count
)
369 set_working_precision(working_precision
);
373 // We want to test the approximation at fixed precision:
374 // either float, double or long double. Begin by getting the
377 boost::math::tools::polynomial
<T
> n
, d
;
378 boost::math::tools::polynomial
<mp_type
> nr
, dr
;
379 nr
= p_remez
->numerator();
380 dr
= p_remez
->denominator();
384 std::vector
<mp_type
> cn1
, cd1
;
385 cn1
= nr
.chebyshev();
386 cd1
= dr
.chebyshev();
387 std::vector
<T
> cn
, cd
;
388 for(unsigned i
= 0; i
< cn1
.size(); ++i
)
390 cn
.push_back(boost::math::tools::real_cast
<T
>(cn1
[i
]));
392 for(unsigned i
= 0; i
< cd1
.size(); ++i
)
394 cd
.push_back(boost::math::tools::real_cast
<T
>(cd1
[i
]));
397 mp_type
max_error(0), max_cheb_error(0);
398 mp_type step
= (b
- a
) / count
;
401 // Do the tests at the zeros:
403 std::cout
<< "Starting tests at " << name
<< " precision...\n";
404 std::cout
<< "Absissa Error (poly) Error (Cheb)\n";
405 for(mp_type x
= a
; x
<= b
; x
+= step
)
407 mp_type true_result
= the_function(x
);
408 //std::cout << true_result << std::endl;
409 T absissa
= boost::math::tools::real_cast
<T
>(x
);
410 mp_type test_result
= convert_to_rr(n
.evaluate(absissa
) / d
.evaluate(absissa
));
411 //std::cout << test_result << std::endl;
412 mp_type cheb_result
= convert_to_rr(boost::math::tools::evaluate_chebyshev(cn
, absissa
) / boost::math::tools::evaluate_chebyshev(cd
, absissa
));
413 //std::cout << cheb_result << std::endl;
414 mp_type err
, cheb_err
;
417 err
= boost::math::tools::relative_error(test_result
, true_result
);
418 cheb_err
= boost::math::tools::relative_error(cheb_result
, true_result
);
422 err
= fabs(test_result
- true_result
);
423 cheb_err
= fabs(cheb_result
- true_result
);
427 if(cheb_err
> max_cheb_error
)
428 max_cheb_error
= cheb_err
;
429 std::cout
<< std::setprecision(6) << std::setw(15) << std::left
<< boost::math::tools::real_cast
<double>(absissa
)
430 << (test_result
< true_result
? "-" : "") << std::setw(20) << std::left
431 << boost::math::tools::real_cast
<double>(err
)
432 << boost::math::tools::real_cast
<double>(cheb_err
) << std::endl
;
434 std::string msg
= "Max Error found at ";
436 msg
+= " precision = ";
437 //msg.append(62 - 17 - msg.size(), ' ');
438 std::cout
<< msg
<< "Poly: " << std::setprecision(6)
439 //<< std::setw(15) << std::left
440 << boost::math::tools::real_cast
<T
>(max_error
)
441 << " Cheb: " << boost::math::tools::real_cast
<T
>(max_cheb_error
) << std::endl
;
445 std::cout
<< "Nothing to test: try converging an approximation first!!!" << std::endl
;
449 void test_n(unsigned n
)
451 do_test_n(mp_type(), "mp_type", n
);
454 void test_float_n(unsigned n
)
456 do_test_n(float(0), "float", n
);
459 void test_double_n(unsigned n
)
461 do_test_n(double(0), "double", n
);
464 void test_long_n(unsigned n
)
466 do_test_n((long double)(0), "long double", n
);
469 void test_float80_n(unsigned n
)
471 do_test_n((boost::multiprecision::cpp_bin_float_double_extended
)(0), "float80", n
);
474 void test_float128_n(unsigned n
)
476 do_test_n((boost::multiprecision::cpp_bin_float_quad
)(0), "float128", n
);
479 void rotate(const char*, const char*)
487 std::cerr
<< "Nothing to rotate" << std::endl
;
491 void rescale(const char*, const char*)
495 p_remez
->rescale(a
, b
);
499 std::cerr
<< "Nothing to rescale" << std::endl
;
503 void graph_poly(const char*, const char*)
506 set_working_precision(working_precision
);
510 // We want to test the approximation at fixed precision:
511 // either float, double or long double. Begin by getting the
514 boost::math::tools::polynomial
<mp_type
> n
, d
;
515 n
= p_remez
->numerator();
516 d
= p_remez
->denominator();
518 mp_type
max_error(0);
519 mp_type step
= (b
- a
) / i
;
521 std::cout
<< "Evaluating Numerator...\n";
523 for(val
= a
; val
<= b
; val
+= step
)
524 std::cout
<< n
.evaluate(val
) << std::endl
;
525 std::cout
<< "Evaluating Denominator...\n";
526 for(val
= a
; val
<= b
; val
+= step
)
527 std::cout
<< d
.evaluate(val
) << std::endl
;
531 std::cout
<< "Nothing to test: try converging an approximation first!!!" << std::endl
;
535 BOOST_AUTO_TEST_CASE( test_main
)
538 real_parser
<long double/*mp_type*/ > const rr_p
;
539 while(std::getline(std::cin
, line
))
541 if(parse(line
.c_str(), str_p("quit"), space_p
).full
)
543 if(false == parse(line
.c_str(),
546 str_p("range")[assign_a(started
, false)] && real_p
[assign_a(a
)] && real_p
[assign_a(b
)]
548 str_p("relative")[assign_a(started
, false)][assign_a(rel_error
, true)]
550 str_p("absolute")[assign_a(started
, false)][assign_a(rel_error
, false)]
552 str_p("pin")[assign_a(started
, false)] && str_p("true")[assign_a(pin
, true)]
554 str_p("pin")[assign_a(started
, false)] && str_p("false")[assign_a(pin
, false)]
556 str_p("pin")[assign_a(started
, false)] && str_p("1")[assign_a(pin
, true)]
558 str_p("pin")[assign_a(started
, false)] && str_p("0")[assign_a(pin
, false)]
560 str_p("pin")[assign_a(started
, false)][assign_a(pin
, true)]
562 str_p("order")[assign_a(started
, false)] && uint_p
[assign_a(orderN
)] && uint_p
[assign_a(orderD
)]
564 str_p("order")[assign_a(started
, false)] && uint_p
[assign_a(orderN
)]
566 str_p("target-precision") && uint_p
[assign_a(target_precision
)]
568 str_p("working-precision")[assign_a(started
, false)] && uint_p
[assign_a(working_precision
)]
570 str_p("variant")[assign_a(started
, false)] && int_p
[assign_a(variant
)]
572 str_p("skew")[assign_a(started
, false)] && int_p
[assign_a(skew
)]
574 str_p("brake") && int_p
[assign_a(brake
)]
576 str_p("step") && int_p
[&step_some
]
580 str_p("poly")[&graph_poly
]
584 str_p("graph") && uint_p
[&do_graph
]
586 str_p("graph")[&graph
]
588 str_p("x-offset") && real_p
[assign_a(x_offset
)]
590 str_p("x-scale") && real_p
[assign_a(x_scale
)]
592 str_p("y-offset") && str_p("auto")[assign_a(auto_offset_y
, true)]
594 str_p("y-offset") && real_p
[assign_a(y_offset
)][assign_a(auto_offset_y
, false)]
596 str_p("test") && str_p("float80") && uint_p
[&test_float80_n
]
598 str_p("test") && str_p("float80")[&test_float80
]
600 str_p("test") && str_p("float128") && uint_p
[&test_float128_n
]
602 str_p("test") && str_p("float128")[&test_float128
]
604 str_p("test") && str_p("float") && uint_p
[&test_float_n
]
606 str_p("test") && str_p("float")[&test_float
]
608 str_p("test") && str_p("double") && uint_p
[&test_double_n
]
610 str_p("test") && str_p("double")[&test_double
]
612 str_p("test") && str_p("long") && uint_p
[&test_long_n
]
614 str_p("test") && str_p("long")[&test_long
]
616 str_p("test") && str_p("all")[&test_all
]
618 str_p("test") && uint_p
[&test_n
]
620 str_p("rotate")[&rotate
]
622 str_p("rescale") && real_p
[assign_a(a
)] && real_p
[assign_a(b
)] && epsilon_p
[&rescale
]
626 std::cout
<< "Unable to parse directive: \"" << line
<< "\"" << std::endl
;
630 std::cout
<< "Variant = " << variant
<< std::endl
;
631 std::cout
<< "range = [" << a
<< "," << b
<< "]" << std::endl
;
632 std::cout
<< "Relative Error = " << rel_error
<< std::endl
;
633 std::cout
<< "Pin to Origin = " << pin
<< std::endl
;
634 std::cout
<< "Order (Num/Denom) = " << orderN
<< "/" << orderD
<< std::endl
;
635 std::cout
<< "Target Precision = " << target_precision
<< std::endl
;
636 std::cout
<< "Working Precision = " << working_precision
<< std::endl
;
637 std::cout
<< "Skew = " << skew
<< std::endl
;
638 std::cout
<< "Brake = " << brake
<< std::endl
;
639 std::cout
<< "X Offset = " << x_offset
<< std::endl
;
640 std::cout
<< "X scale = " << x_scale
<< std::endl
;
641 std::cout
<< "Y Offset = ";
643 std::cout
<< "Auto (";
644 std::cout
<< y_offset
;
647 std::cout
<< std::endl
;