mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 04:22:09 +00:00
Reinch's modification to Clenshaw recurrence (#339)
* Reinch's modification to Clenshaw recurrence. [CI SKIP] * Convert Chebyshev tests to math_unit_test.hpp * Performance of translated Chebyshev Clenshaw recurrence. [CI SKIP] * Prepare to use modified Clenshaw recurrence in Chebyshev transform. * Remove unused headers from Chebyshev transform test [CI SKIP] * Update Chebyshev transform tests to use math_unit_test.hpp
This commit is contained in:
Nick
GitHub
109
doc/fp_utilities/ulps_plots.qbk
Normal file
109
doc/fp_utilities/ulps_plots.qbk
Normal file
@@ -0,0 +1,109 @@
|
||||
[/
|
||||
Copyright Nick Thompson, John Maddock, 2020
|
||||
Distributed under the Boost Software License, Version 1.0.
|
||||
(See accompanying file LICENSE_1_0.txt or copy at
|
||||
http://www.boost.org/LICENSE_1_0.txt).
|
||||
]
|
||||
|
||||
[section:ulps_plots ULPs Plots]
|
||||
|
||||
[h4 Synopsis]
|
||||
|
||||
#include <boost/math/tools/ulps_plot.hpp>
|
||||
|
||||
namespace boost::math::tools {
|
||||
|
||||
template<class F, typename PreciseReal, typename CoarseReal>
|
||||
class ulps_plot {
|
||||
public:
|
||||
ulps_plot(F hi_acc_impl, CoarseReal a, CoarseReal b,
|
||||
size_t samples = 10000, bool perturb_abscissas = false, int random_seed = -1);
|
||||
|
||||
ulps_plot& clip(PreciseReal clip);
|
||||
|
||||
ulps_plot& width(int width);
|
||||
|
||||
ulps_plot& envelope_color(std::string const & color);
|
||||
|
||||
ulps_plot& title(std::string const & title);
|
||||
|
||||
ulps_plot& background_color(std::string const & background_color);
|
||||
|
||||
ulps_plot& font_color(std::string const & font_color);
|
||||
|
||||
ulps_plot& ulp_envelope(bool write_ulp);
|
||||
|
||||
ulps_plot& horizontal_lines(int horizontal_lines);
|
||||
|
||||
ulps_plot& vertical_lines(int vertical_lines);
|
||||
|
||||
template<class G>
|
||||
ulps_plot& add_fn(G g, std::string const & color = "steelblue");
|
||||
|
||||
friend std::ostream& operator<<(std::ostream& fs, ulps_plot const & plot)
|
||||
|
||||
void write(std::string const & filename) const;
|
||||
|
||||
};
|
||||
|
||||
} // namespaces
|
||||
|
||||
|
||||
An ULPs plot measures the accuracy of an implementation of a function implemented in floating point arithmetic.
|
||||
ULP stands for /unit in the last place/; i.e., we create a unit system that normalizes the distance between one representable and the next greater representable to 1.
|
||||
|
||||
So for example, the ULP distance between `1.0` and `1.0 + std::numeric_limits<double>::epsilon()` is 1, the ulp distance between `1.0` and `1.0 + 2*std::numeric_limits<double>::epsilon()` is 2.
|
||||
A key concept of this measurement is /semi-scale invariance/: The ULP distance between `2.0` and `2.0 + 2*std::numeric_limits<double>::epsilon()` is 1, not 2.
|
||||
|
||||
An ULPs plot gives us an idea of how efficiently we are using the floating point number system we have chosen.
|
||||
To compute a ULP plot, we need a higher precision implementation; for example we can test a 32 bit floating point implementation against an 80-bit long double implementation.
|
||||
The higher precision value is computed in the template type `PreciseReal`, and the lower precision value is computed in type `CoarseReal`.
|
||||
For simplicity, we will refer to the result of the higher precision implementation the exact value, which we will denote by /y/,
|
||||
and the value computed in lower precision we will denote by ŷ.
|
||||
The ULP distance between /y/ and ŷ is then
|
||||
|
||||
[$../graphs/ulp_distance.svg]
|
||||
|
||||
where μ is the unit roundoff of the `CoarseReal` type.
|
||||
|
||||
If every sample in an ULPs plot is bounded by ±½, then we have good evidence that we have used our floating point number as efficiently as we can possibly expect:
|
||||
We are getting the correctly rounded result.
|
||||
However, the converse is not true: If there exists samples that are well above the unit roundoff, we do not have evidence that the function is poorly implemented.
|
||||
Considering the /y/=0 case makes this obvious, but another more subtle issue is at play.
|
||||
Consider that when we evaluate a function /f/, it is rare that we want to evaluate /f/ at a representable x̂,
|
||||
instead we have an infinite precision value /x/ which must be rounded into x̂ using the floating point rounding model x̂ = x(1+ε), where |ε|<μ.
|
||||
|
||||
We can use a first order Taylor series to determine the best accuracy we can hope for under the assumption that the only source of error is rounding the abscissa to the nearest to nearest representable:
|
||||
|
||||
[$../graphs/ulps_and_condition_number.svg]
|
||||
|
||||
Hence the condition number of function evaluation provides a natural scale for evaluating the quality of an implementation of a special function.
|
||||
Since, in addition, we cannot expect better than half-ULP accuracy, we can create a /ulp envelope/ by taking the maximum of 0.5 and half the condition number of function evaluation at each point.
|
||||
|
||||
[$../graphs/airy_ai_float.svg]
|
||||
|
||||
The graph shows the ULP accuracy of Boost's implementation of the Airy Ai function.
|
||||
We see only a few evaluations outside the ULP envelope.
|
||||
Note how the condition number of function evaluation becomes unbounded at the roots.
|
||||
|
||||
|
||||
A minimal working example which reproduces the plot above is given [@../../example/airy_ulps_plot.cpp here].
|
||||
|
||||
Note that you can compare two different implementations of a single function by calling `add_fn` multiple times.
|
||||
This will give you a graphical answer to which one is more accurate.
|
||||
For example:
|
||||
|
||||
auto plot = ulps_plot<decltype(f), long double, float>(f, 1.0f, 2.0f);
|
||||
plot.add_fn(impl1, "steelblue");
|
||||
plot.add_fn(impl2, "orange");
|
||||
|
||||
|
||||
Note that the graph is written as an SVG. We use `firefox foo.svg` or `open -a Firefox foo.svg` to display these files.
|
||||
|
||||
[heading References]
|
||||
|
||||
* Cleve Moler. ['Ulps Plots Reveal Math Function Accuracy], https://blogs.mathworks.com/cleve/2017/01/23/ulps-plots-reveal-math-function-accurary/
|
||||
|
||||
|
||||
|
||||
[endsect]
|
||||
9917
doc/graphs/airy_ai_float.svg
Normal file
9917
doc/graphs/airy_ai_float.svg
Normal file
File diff suppressed because one or more lines are too long
|
After Width: | Height: | Size: 758 KiB |
2
doc/graphs/chebyshev_convention_1.svg
Normal file
2
doc/graphs/chebyshev_convention_1.svg
Normal file
File diff suppressed because one or more lines are too long
|
After Width: | Height: | Size: 8.9 KiB |
2
doc/graphs/chebyshev_convention_2.svg
Normal file
2
doc/graphs/chebyshev_convention_2.svg
Normal file
File diff suppressed because one or more lines are too long
|
After Width: | Height: | Size: 10 KiB |
2
doc/graphs/ulp_distance.svg
Normal file
2
doc/graphs/ulp_distance.svg
Normal file
@@ -0,0 +1,2 @@
|
||||
<?xml version="1.0" encoding="UTF-8" standalone="no" ?>
|
||||
<svg xmlns="http://www.w3.org/2000/svg" width="49.016px" height="43.440px" viewBox="0 -1441 2708 2400" xmlns:xlink="http://www.w3.org/1999/xlink"><defs><path id="MJX-24-TEX-I-79" d="M21 287Q21 301 36 335T84 406T158 442Q199 442 224 419T250 355Q248 336 247 334Q247 331 231 288T198 191T182 105Q182 62 196 45T238 27Q261 27 281 38T312 61T339 94Q339 95 344 114T358 173T377 247Q415 397 419 404Q432 431 462 431Q475 431 483 424T494 412T496 403Q496 390 447 193T391 -23Q363 -106 294 -155T156 -205Q111 -205 77 -183T43 -117Q43 -95 50 -80T69 -58T89 -48T106 -45Q150 -45 150 -87Q150 -107 138 -122T115 -142T102 -147L99 -148Q101 -153 118 -160T152 -167H160Q177 -167 186 -165Q219 -156 247 -127T290 -65T313 -9T321 21L315 17Q309 13 296 6T270 -6Q250 -11 231 -11Q185 -11 150 11T104 82Q103 89 103 113Q103 170 138 262T173 379Q173 380 173 381Q173 390 173 393T169 400T158 404H154Q131 404 112 385T82 344T65 302T57 280Q55 278 41 278H27Q21 284 21 287Z"></path><path id="MJX-24-TEX-N-5E" d="M112 560L249 694L257 686Q387 562 387 560L361 531Q359 532 303 581L250 627L195 580Q182 569 169 557T148 538L140 532Q138 530 125 546L112 560Z"></path><path id="MJX-24-TEX-N-2212" d="M84 237T84 250T98 270H679Q694 262 694 250T679 230H98Q84 237 84 250Z"></path><path id="MJX-24-TEX-N-32" d="M109 429Q82 429 66 447T50 491Q50 562 103 614T235 666Q326 666 387 610T449 465Q449 422 429 383T381 315T301 241Q265 210 201 149L142 93L218 92Q375 92 385 97Q392 99 409 186V189H449V186Q448 183 436 95T421 3V0H50V19V31Q50 38 56 46T86 81Q115 113 136 137Q145 147 170 174T204 211T233 244T261 278T284 308T305 340T320 369T333 401T340 431T343 464Q343 527 309 573T212 619Q179 619 154 602T119 569T109 550Q109 549 114 549Q132 549 151 535T170 489Q170 464 154 447T109 429Z"></path><path id="MJX-24-TEX-I-3BC" d="M58 -216Q44 -216 34 -208T23 -186Q23 -176 96 116T173 414Q186 442 219 442Q231 441 239 435T249 423T251 413Q251 401 220 279T187 142Q185 131 185 107V99Q185 26 252 26Q261 26 270 27T287 31T302 38T315 45T327 55T338 65T348 77T356 88T365 100L372 110L408 253Q444 395 448 404Q461 431 491 431Q504 431 512 424T523 412T525 402L449 84Q448 79 448 68Q448 43 455 35T476 26Q485 27 496 35Q517 55 537 131Q543 151 547 152Q549 153 557 153H561Q580 153 580 144Q580 138 575 117T555 63T523 13Q510 0 491 -8Q483 -10 467 -10Q446 -10 429 -4T402 11T385 29T376 44T374 51L368 45Q362 39 350 30T324 12T288 -4T246 -11Q199 -11 153 12L129 -85Q108 -167 104 -180T92 -202Q76 -216 58 -216Z"></path><path id="MJX-24-TEX-N-7C" d="M139 -249H137Q125 -249 119 -235V251L120 737Q130 750 139 750Q152 750 159 735V-235Q151 -249 141 -249H139Z"></path></defs><g stroke="currentColor" fill="currentColor" stroke-width="0" transform="matrix(1 0 0 -1 0 0)"><g data-mml-node="math"><g data-mml-node="mfrac"><g data-mml-node="mrow" transform="translate(220, 676)"><g data-mml-node="TeXAtom"><g data-mml-node="mover"><g data-mml-node="mi" transform="translate(5, 0)"><use xlink:href="#MJX-24-TEX-I-79"></use></g><g data-mml-node="mo" transform="translate(55.6, -29)"><use xlink:href="#MJX-24-TEX-N-5E"></use></g></g></g><g data-mml-node="mo" transform="translate(777.8, 0)"><use xlink:href="#MJX-24-TEX-N-2212"></use></g><g data-mml-node="mi" transform="translate(1778, 0)"><use xlink:href="#MJX-24-TEX-I-79"></use></g></g><g data-mml-node="mrow" transform="translate(279.5, -709.5)"><g data-mml-node="mn"><use xlink:href="#MJX-24-TEX-N-32"></use></g><g data-mml-node="mi" transform="translate(500, 0)"><use xlink:href="#MJX-24-TEX-I-3BC"></use></g><g data-mml-node="TeXAtom" transform="translate(1103, 0)"><g data-mml-node="mo"><use xlink:href="#MJX-24-TEX-N-7C"></use></g></g><g data-mml-node="mi" transform="translate(1381, 0)"><use xlink:href="#MJX-24-TEX-I-79"></use></g><g data-mml-node="mo" transform="translate(1871, 0)"><use xlink:href="#MJX-24-TEX-N-7C"></use></g></g><rect width="2468" height="60" x="120" y="220"></rect></g></g></g></svg>
|
||||
|
After Width: | Height: | Size: 3.8 KiB |
2
doc/graphs/ulps_and_condition_number.svg
Normal file
2
doc/graphs/ulps_and_condition_number.svg
Normal file
File diff suppressed because one or more lines are too long
|
After Width: | Height: | Size: 19 KiB |
@@ -574,6 +574,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22.
|
||||
[include fp_utilities/float_next.qbk]
|
||||
[include fp_utilities/float_comparison.qbk]
|
||||
[include fp_utilities/condition_numbers.qbk]
|
||||
[include fp_utilities/ulps_plots.qbk]
|
||||
[endmathpart]
|
||||
|
||||
[mathpart cstdfloat..Specified-width floating-point typedefs]
|
||||
|
||||
@@ -34,7 +34,10 @@
|
||||
``__sf_result`` chebyshev_t_prime(unsigned n, Real const & x);
|
||||
|
||||
template<class Real1, class Real2>
|
||||
``__sf_result`` chebyshev_clenshaw_recurrence(const Real* const c, size_t length, Real2 x);
|
||||
``__sf_result`` chebyshev_clenshaw_recurrence(const Real1* const c, size_t length, Real2 x);
|
||||
|
||||
template<typename Real>
|
||||
Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, Real a, Real b, Real x);
|
||||
|
||||
}} // namespaces
|
||||
|
||||
@@ -78,20 +81,32 @@ and is implemented in boost as
|
||||
|
||||
double x = 0.5;
|
||||
std::vector<double> c{14.2, -13.7, 82.3, 96};
|
||||
double f = chebyshev_clenshaw_recurrence(c.data(), c.size(), Real x);
|
||||
double f = chebyshev_clenshaw_recurrence(c.data(), c.size(), x);
|
||||
|
||||
|
||||
N.B.: There is factor of /2/ difference in our definition of the first coefficient in the Chebyshev series from Clenshaw's original work.
|
||||
This is because two traditions exist in notation for the Chebyshev series expansion,
|
||||
|
||||
[:/f/(/x/) \u2248 \u2211[sub n=0][super N-1] /a/[sub n]/T/[sub n](/x/)]
|
||||
[$../graphs/chebyshev_convention_1.svg]
|
||||
|
||||
|
||||
and
|
||||
|
||||
[:/f/(/x/) \u2248 /c/[sub 0]/2 + \u2211[sub n=1][super N-1] /c/[sub n]/T/[sub n](/x/)]
|
||||
[$../graphs/chebyshev_convention_2.svg]
|
||||
|
||||
['*boost math always uses the second convention, with the factor of 1/2 on the first coefficient.*]
|
||||
|
||||
Another common use case is when the polynomial must be evaluated on some interval \[/a/, /b/\].
|
||||
The translation to the interval \[-1, 1\] causes a few accuracy problems and also gives us some opportunities.
|
||||
For this case, we use [@https://doi.org/10.1093/imamat/20.3.379 Reinch's modification] to the Clenshaw recurrence, which is also discussed [@https://archive.siam.org/books/ot99/OT99SampleChapter.pdf here].
|
||||
The usage is as follows:
|
||||
|
||||
double x = 9;
|
||||
double a = 7;
|
||||
double b = 12;
|
||||
std::vector<double> c{14.2, -13.7, 82.3, 96};
|
||||
double f = chebyshev_clenshaw_recurrence(c.data(), c.size(), a, b, x);
|
||||
|
||||
Chebyshev polynomials of the second kind can be evaluated via `chebyshev_u`:
|
||||
|
||||
double x = -0.23;
|
||||
@@ -151,7 +166,7 @@ The notion of "very close" can be made rigorous; see Trefethen's "Approximation
|
||||
|
||||
The Chebyshev transform works by creating a vector of values by evaluating the input function at the Chebyshev points, and then performing a discrete cosine transform on the resulting vector.
|
||||
In order to do this efficiently, we have used [@http://www.fftw.org/ FFTW3].
|
||||
So to compile, you must have `FFTW3` installed, and link with `-lfftw3` for double precision, `-lfftw3f` for float precision, `-lfftw3l` for long double precision, and -lfftwq for quad (`__float128`) precision.
|
||||
So to compile, you must have `FFTW3` installed, and link with `-lfftw3` for double precision, `-lfftw3f` for float precision, `-lfftw3l` for long double precision, and `-lfftw3q` for quad (`__float128`) precision.
|
||||
After the coefficients of the Chebyshev series are known, the routine goes back through them and filters out all the coefficients whose absolute ratio to the largest coefficient are less than the tolerance requested in the constructor.
|
||||
|
||||
[endsect] [/section:chebyshev Chebyshev Polynomials]
|
||||
|
||||
@@ -73,6 +73,7 @@ test-suite examples :
|
||||
[ run hyperexponential_more_snips.cpp ]
|
||||
[ run inverse_chi_squared_example.cpp ]
|
||||
[ run legendre_stieltjes_example.cpp : : : [ requires cxx11_auto_declarations cxx11_defaulted_functions cxx11_lambdas ] ]
|
||||
[ run airy_ulps_plot.cpp : : : [ requires cxx17_std_apply cxx17_if_constexpr ] ]
|
||||
#[ # run inverse_chi_squared_find_df_example.cpp ]
|
||||
#[ run lambert_w_basic_example.cpp ]
|
||||
[ run lambert_w_basic_example.cpp : : : [ requires cxx11_numeric_limits ] ]
|
||||
|
||||
60
example/airy_ulps_plot.cpp
Normal file
60
example/airy_ulps_plot.cpp
Normal file
@@ -0,0 +1,60 @@
|
||||
// (C) Copyright Nick Thompson 2020.
|
||||
// Use, modification and distribution are subject to the
|
||||
// Boost Software License, Version 1.0. (See accompanying file
|
||||
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||||
#include <iostream>
|
||||
#include <boost/math/tools/ulps_plot.hpp>
|
||||
#include <boost/core/demangle.hpp>
|
||||
#include <boost/math/special_functions/airy.hpp>
|
||||
|
||||
using boost::math::tools::ulps_plot;
|
||||
|
||||
int main() {
|
||||
using PreciseReal = long double;
|
||||
using CoarseReal = float;
|
||||
|
||||
typedef boost::math::policies::policy<
|
||||
boost::math::policies::promote_float<false>,
|
||||
boost::math::policies::promote_double<false> >
|
||||
no_promote_policy;
|
||||
|
||||
auto ai_coarse = [](CoarseReal x) {
|
||||
return boost::math::airy_ai<CoarseReal>(x, no_promote_policy());
|
||||
};
|
||||
auto ai_precise = [](PreciseReal x) {
|
||||
return boost::math::airy_ai<PreciseReal>(x, no_promote_policy());
|
||||
};
|
||||
|
||||
std::string filename = "airy_ai_" + boost::core::demangle(typeid(CoarseReal).name()) + ".svg";
|
||||
int samples = 10000;
|
||||
// How many pixels wide do you want your .svg?
|
||||
int width = 700;
|
||||
// Near a root, we have unbounded relative error. So for functions with roots, we define an ULP clip:
|
||||
PreciseReal clip = 2.5;
|
||||
// Should we perturb the abscissas? i.e., should we compute the high precision function f at x,
|
||||
// and the low precision function at the nearest representable x̂ to x?
|
||||
// Or should we compute both the high precision and low precision function at a low precision representable x̂?
|
||||
bool perturb_abscissas = false;
|
||||
auto plot = ulps_plot<decltype(ai_precise), PreciseReal, CoarseReal>(ai_precise, CoarseReal(-3), CoarseReal(3), samples, perturb_abscissas);
|
||||
// Note the argument chaining:
|
||||
plot.clip(clip).width(width);
|
||||
// Sometimes it's useful to set a title, but in many cases it's more useful to just use a caption.
|
||||
//std::string title = "Airy Ai ULP plot at " + boost::core::demangle(typeid(CoarseReal).name()) + " precision";
|
||||
//plot.title(title);
|
||||
plot.vertical_lines(6);
|
||||
plot.add_fn(ai_coarse);
|
||||
// You can write the plot to a stream:
|
||||
//std::cout << plot;
|
||||
// Or to a file:
|
||||
plot.write(filename);
|
||||
|
||||
// Don't like the default dark theme?
|
||||
plot.background_color("white").font_color("black");
|
||||
filename = "airy_ai_" + boost::core::demangle(typeid(CoarseReal).name()) + "_white.svg";
|
||||
plot.write(filename);
|
||||
|
||||
// Don't like the envelope?
|
||||
plot.ulp_envelope(false);
|
||||
filename = "airy_ai_" + boost::core::demangle(typeid(CoarseReal).name()) + "_white_no_envelope.svg";
|
||||
plot.write(filename);
|
||||
}
|
||||
@@ -11,8 +11,8 @@
|
||||
|
||||
#include <boost/multiprecision/float128.hpp>
|
||||
#include <boost/math/special_functions/daubechies_scaling.hpp>
|
||||
#include <boost/math/tools/ulp_plot.hpp>
|
||||
#include <quicksvg/graph_fn.hpp>
|
||||
#include <quicksvg/ulp_plot.hpp>
|
||||
|
||||
|
||||
using boost::multiprecision::float128;
|
||||
@@ -135,11 +135,16 @@ void do_ulp(int coarse_refinements, PhiPrecise phi_precise)
|
||||
title = "";
|
||||
|
||||
std::string filename = "daubechies_" + std::to_string(p) + "_" + boost::core::demangle(typeid(CoarseReal).name()) + "_" + std::to_string(coarse_refinements) + "_refinements.svg";
|
||||
int samples = 20000;
|
||||
int samples = 1000000;
|
||||
int clip = 20;
|
||||
int horizontal_lines = 8;
|
||||
int vertical_lines = 2*p - 1;
|
||||
quicksvg::ulp_plot<decltype(phi_coarse), CoarseReal, decltype(phi_precise), PreciseReal>(phi_coarse, phi_precise, CoarseReal(0), phi_coarse.support().second, title, filename, samples, GRAPH_WIDTH, clip, horizontal_lines, vertical_lines);
|
||||
auto [a, b] = phi_coarse.support();
|
||||
auto plot = boost::math::tools::ulp_plot<decltype(phi_precise), PreciseReal, CoarseReal>(phi_precise, a, b, samples);
|
||||
plot.set_clip(clip);
|
||||
plot.set_width(GRAPH_WIDTH);
|
||||
plot.add_fn(phi_coarse);
|
||||
plot.write(filename, true, title, horizontal_lines, vertical_lines);
|
||||
}
|
||||
|
||||
|
||||
@@ -151,9 +156,9 @@ int main()
|
||||
boost::hana::for_each(std::make_index_sequence<18>(), [&](auto i){ plot_convergence<double, i+2>(); });
|
||||
|
||||
using PreciseReal = float128;
|
||||
using CoarseReal = double;
|
||||
using CoarseReal = float;
|
||||
int precise_refinements = 22;
|
||||
constexpr const int p = 8;
|
||||
constexpr const int p = 14;
|
||||
std::cout << "Computing precise scaling function in " << boost::core::demangle(typeid(PreciseReal).name()) << " precision.\n";
|
||||
auto phi_precise = boost::math::daubechies_scaling<PreciseReal, p>(precise_refinements);
|
||||
std::cout << "Beginning comparison with functions computed in " << boost::core::demangle(typeid(CoarseReal).name()) << " precision.\n";
|
||||
|
||||
@@ -165,5 +165,100 @@ inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, co
|
||||
}
|
||||
|
||||
|
||||
|
||||
namespace detail {
|
||||
template<class Real>
|
||||
inline Real unchecked_chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real & a, const Real & b, const Real& x)
|
||||
{
|
||||
Real t;
|
||||
Real u;
|
||||
// This cutoff is not super well defined, but it's a good estimate.
|
||||
// See "An Error Analysis of the Modified Clenshaw Method for Evaluating Chebyshev and Fourier Series"
|
||||
// J. OLIVER, IMA Journal of Applied Mathematics, Volume 20, Issue 3, November 1977, Pages 379–391
|
||||
// https://doi.org/10.1093/imamat/20.3.379
|
||||
const Real cutoff = 0.6;
|
||||
if (x - a < b - x)
|
||||
{
|
||||
u = 2*(x-a)/(b-a);
|
||||
t = u - 1;
|
||||
if (t > -cutoff)
|
||||
{
|
||||
Real b2 = 0;
|
||||
Real b1 = c[length -1];
|
||||
for(size_t j = length - 2; j >= 1; --j)
|
||||
{
|
||||
Real tmp = 2*t*b1 - b2 + c[j];
|
||||
b2 = b1;
|
||||
b1 = tmp;
|
||||
}
|
||||
return t*b1 - b2 + c[0]/2;
|
||||
}
|
||||
else
|
||||
{
|
||||
Real b = c[length -1];
|
||||
Real d = b;
|
||||
Real b2 = 0;
|
||||
for (size_t r = length - 2; r >= 1; --r)
|
||||
{
|
||||
d = 2*u*b - d + c[r];
|
||||
b2 = b;
|
||||
b = d - b;
|
||||
}
|
||||
return t*b - b2 + c[0]/2;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
u = -2*(b-x)/(b-a);
|
||||
t = u + 1;
|
||||
if (t < cutoff)
|
||||
{
|
||||
Real b2 = 0;
|
||||
Real b1 = c[length -1];
|
||||
for(size_t j = length - 2; j >= 1; --j)
|
||||
{
|
||||
Real tmp = 2*t*b1 - b2 + c[j];
|
||||
b2 = b1;
|
||||
b1 = tmp;
|
||||
}
|
||||
return t*b1 - b2 + c[0]/2;
|
||||
}
|
||||
else
|
||||
{
|
||||
Real b = c[length -1];
|
||||
Real d = b;
|
||||
Real b2 = 0;
|
||||
for (size_t r = length - 2; r >= 1; --r)
|
||||
{
|
||||
d = 2*u*b + d + c[r];
|
||||
b2 = b;
|
||||
b = d + b;
|
||||
}
|
||||
return t*b - b2 + c[0]/2;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
|
||||
template<class Real>
|
||||
inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real & a, const Real & b, const Real& x)
|
||||
{
|
||||
if (x < a || x > b)
|
||||
{
|
||||
throw std::domain_error("x in [a, b] is required.");
|
||||
}
|
||||
if (length < 2)
|
||||
{
|
||||
if (length == 0)
|
||||
{
|
||||
return 0;
|
||||
}
|
||||
return c[0]/2;
|
||||
}
|
||||
return detail::unchecked_chebyshev_clenshaw_recurrence(c, length, a, b, x);
|
||||
}
|
||||
|
||||
|
||||
}}
|
||||
#endif
|
||||
|
||||
@@ -116,7 +116,7 @@ public:
|
||||
template<class F>
|
||||
chebyshev_transform(const F& f, Real a, Real b,
|
||||
Real tol = 500 * std::numeric_limits<Real>::epsilon(),
|
||||
size_t max_refinements = 15) : m_a(a), m_b(b)
|
||||
size_t max_refinements = 16) : m_a(a), m_b(b)
|
||||
{
|
||||
if (a >= b)
|
||||
{
|
||||
@@ -174,16 +174,9 @@ public:
|
||||
}
|
||||
}
|
||||
|
||||
Real operator()(Real x) const
|
||||
inline Real operator()(Real x) const
|
||||
{
|
||||
using boost::math::constants::half;
|
||||
if (x > m_b || x < m_a)
|
||||
{
|
||||
throw std::domain_error("x not in [a, b]\n");
|
||||
}
|
||||
|
||||
Real z = (2*x - m_a - m_b)/(m_b - m_a);
|
||||
return chebyshev_clenshaw_recurrence(m_coeffs.data(), m_coeffs.size(), z);
|
||||
return chebyshev_clenshaw_recurrence(m_coeffs.data(), m_coeffs.size(), m_a, m_b, x);
|
||||
}
|
||||
|
||||
// Integral over entire domain [a, b]
|
||||
|
||||
541
include/boost/math/tools/ulps_plot.hpp
Normal file
541
include/boost/math/tools/ulps_plot.hpp
Normal file
@@ -0,0 +1,541 @@
|
||||
// (C) Copyright Nick Thompson 2020.
|
||||
// Use, modification and distribution are subject to the
|
||||
// Boost Software License, Version 1.0. (See accompanying file
|
||||
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||||
#ifndef BOOST_MATH_TOOLS_ULP_PLOT_HPP
|
||||
#define BOOST_MATH_TOOLS_ULP_PLOT_HPP
|
||||
#include <algorithm>
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
#include <cassert>
|
||||
#include <vector>
|
||||
#include <utility>
|
||||
#include <fstream>
|
||||
#include <string>
|
||||
#include <list>
|
||||
#include <random>
|
||||
#include <stdexcept>
|
||||
#include <boost/math/tools/condition_numbers.hpp>
|
||||
#include <boost/random/uniform_real_distribution.hpp>
|
||||
#include <boost/algorithm/string/predicate.hpp>
|
||||
|
||||
|
||||
// Design of this function comes from:
|
||||
// https://blogs.mathworks.com/cleve/2017/01/23/ulps-plots-reveal-math-function-accurary/
|
||||
|
||||
// The envelope is the maximum of 1/2 and half the condition number of function evaluation.
|
||||
|
||||
namespace boost::math::tools {
|
||||
|
||||
namespace detail {
|
||||
template<class F1, class F2, class Real>
|
||||
void write_gridlines(std::ostream& fs, int horizontal_lines, int vertical_lines,
|
||||
F1 x_scale, F2 y_scale, Real min_x, Real max_x, Real min_y, Real max_y,
|
||||
int graph_width, int graph_height, int margin_left, std::string const & font_color)
|
||||
{
|
||||
// Make a grid:
|
||||
for (int i = 1; i <= horizontal_lines; ++i) {
|
||||
Real y_cord_dataspace = min_y + ((max_y - min_y)*i)/horizontal_lines;
|
||||
auto y = y_scale(y_cord_dataspace);
|
||||
fs << "<line x1='0' y1='" << y << "' x2='" << graph_width
|
||||
<< "' y2='" << y
|
||||
<< "' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />\n";
|
||||
|
||||
fs << "<text x='" << -margin_left/4 + 5 << "' y='" << y - 3
|
||||
<< "' font-family='times' font-size='10' fill='" << font_color << "' transform='rotate(-90 "
|
||||
<< -margin_left/4 + 8 << " " << y + 5 << ")'>"
|
||||
<< std::setprecision(4) << y_cord_dataspace << "</text>\n";
|
||||
}
|
||||
|
||||
for (int i = 1; i <= vertical_lines; ++i) {
|
||||
Real x_cord_dataspace = min_x + ((max_x - min_x)*i)/vertical_lines;
|
||||
Real x = x_scale(x_cord_dataspace);
|
||||
fs << "<line x1='" << x << "' y1='0' x2='" << x
|
||||
<< "' y2='" << graph_height
|
||||
<< "' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />\n";
|
||||
|
||||
fs << "<text x='" << x - 10 << "' y='" << graph_height + 10
|
||||
<< "' font-family='times' font-size='10' fill='" << font_color << "'>"
|
||||
<< std::setprecision(4) << x_cord_dataspace << "</text>\n";
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<class F, typename PreciseReal, typename CoarseReal>
|
||||
class ulps_plot {
|
||||
public:
|
||||
ulps_plot(F hi_acc_impl, CoarseReal a, CoarseReal b,
|
||||
size_t samples = 10000, bool perturb_abscissas = false, int random_seed = -1);
|
||||
|
||||
ulps_plot& clip(PreciseReal clip);
|
||||
|
||||
ulps_plot& width(int width);
|
||||
|
||||
ulps_plot& envelope_color(std::string const & color);
|
||||
|
||||
ulps_plot& title(std::string const & title);
|
||||
|
||||
ulps_plot& background_color(std::string const & background_color);
|
||||
|
||||
ulps_plot& font_color(std::string const & font_color);
|
||||
|
||||
ulps_plot& ulp_envelope(bool write_ulp);
|
||||
|
||||
template<class G>
|
||||
ulps_plot& add_fn(G g, std::string const & color = "steelblue");
|
||||
|
||||
ulps_plot& horizontal_lines(int horizontal_lines);
|
||||
|
||||
ulps_plot& vertical_lines(int vertical_lines);
|
||||
|
||||
void write(std::string const & filename) const;
|
||||
|
||||
friend std::ostream& operator<<(std::ostream& fs, ulps_plot const & plot)
|
||||
{
|
||||
using std::abs;
|
||||
using std::floor;
|
||||
using std::isnan;
|
||||
if (plot.ulp_list_.size() == 0)
|
||||
{
|
||||
throw std::domain_error("No functions added for comparison.");
|
||||
}
|
||||
if (plot.width_ <= 1)
|
||||
{
|
||||
throw std::domain_error("Width = " + std::to_string(plot.width_) + ", which is too small.");
|
||||
}
|
||||
|
||||
PreciseReal worst_ulp_distance = 0;
|
||||
PreciseReal min_y = std::numeric_limits<PreciseReal>::max();
|
||||
PreciseReal max_y = std::numeric_limits<PreciseReal>::lowest();
|
||||
for (auto const & ulp_vec : plot.ulp_list_)
|
||||
{
|
||||
for (auto const & ulp : ulp_vec)
|
||||
{
|
||||
if (abs(ulp) > worst_ulp_distance)
|
||||
{
|
||||
worst_ulp_distance = abs(ulp);
|
||||
}
|
||||
if (ulp < min_y)
|
||||
{
|
||||
min_y = ulp;
|
||||
}
|
||||
if (ulp > max_y)
|
||||
{
|
||||
max_y = ulp;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// half-ulp accuracy is the best that can be expected; sometimes we can get less, but barely less.
|
||||
// then the axes don't show up; painful!
|
||||
if (max_y < 0.5) {
|
||||
max_y = 0.5;
|
||||
}
|
||||
if (min_y > -0.5) {
|
||||
min_y = -0.5;
|
||||
}
|
||||
|
||||
if (plot.clip_ > 0)
|
||||
{
|
||||
if (max_y > plot.clip_)
|
||||
{
|
||||
max_y = plot.clip_;
|
||||
}
|
||||
if (min_y < -plot.clip_)
|
||||
{
|
||||
min_y = -plot.clip_;
|
||||
}
|
||||
}
|
||||
|
||||
int height = floor(double(plot.width_)/1.61803);
|
||||
int margin_top = 40;
|
||||
int margin_left = 25;
|
||||
if (plot.title_.size() == 0)
|
||||
{
|
||||
margin_top = 10;
|
||||
margin_left = 15;
|
||||
}
|
||||
int margin_bottom = 20;
|
||||
int margin_right = 20;
|
||||
int graph_height = height - margin_bottom - margin_top;
|
||||
int graph_width = plot.width_ - margin_left - margin_right;
|
||||
|
||||
// Maps [a,b] to [0, graph_width]
|
||||
auto x_scale = [&](CoarseReal x)->CoarseReal
|
||||
{
|
||||
return ((x-plot.a_)/(plot.b_ - plot.a_))*static_cast<CoarseReal>(graph_width);
|
||||
};
|
||||
|
||||
auto y_scale = [&](PreciseReal y)->PreciseReal
|
||||
{
|
||||
return ((max_y - y)/(max_y - min_y) )*static_cast<PreciseReal>(graph_height);
|
||||
};
|
||||
|
||||
fs << "<?xml version=\"1.0\" encoding='UTF-8' ?>\n"
|
||||
<< "<svg xmlns='http://www.w3.org/2000/svg' width='"
|
||||
<< plot.width_ << "' height='"
|
||||
<< height << "'>\n"
|
||||
<< "<style>\nsvg { background-color:" << plot.background_color_ << "; }\n"
|
||||
<< "</style>\n";
|
||||
if (plot.title_.size() > 0)
|
||||
{
|
||||
fs << "<text x='" << floor(plot.width_/2)
|
||||
<< "' y='" << floor(margin_top/2)
|
||||
<< "' font-family='Palatino' font-size='25' fill='"
|
||||
<< plot.font_color_ << "' alignment-baseline='middle' text-anchor='middle'>"
|
||||
<< plot.title_
|
||||
<< "</text>\n";
|
||||
}
|
||||
|
||||
// Construct SVG group to simplify the calculations slightly:
|
||||
fs << "<g transform='translate(" << margin_left << ", " << margin_top << ")'>\n";
|
||||
// y-axis:
|
||||
fs << "<line x1='0' y1='0' x2='0' y2='" << graph_height
|
||||
<< "' stroke='gray' stroke-width='1'/>\n";
|
||||
PreciseReal x_axis_loc = y_scale(static_cast<PreciseReal>(0));
|
||||
fs << "<line x1='0' y1='" << x_axis_loc
|
||||
<< "' x2='" << graph_width << "' y2='" << x_axis_loc
|
||||
<< "' stroke='gray' stroke-width='1'/>\n";
|
||||
|
||||
if (worst_ulp_distance > 3)
|
||||
{
|
||||
detail::write_gridlines(fs, plot.horizontal_lines_, plot.vertical_lines_, x_scale, y_scale, plot.a_, plot.b_,
|
||||
static_cast<CoarseReal>(min_y), static_cast<CoarseReal>(max_y), graph_width, graph_height, margin_left, plot.font_color_);
|
||||
}
|
||||
else
|
||||
{
|
||||
std::vector<double> ys{-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0};
|
||||
for (size_t i = 0; i < ys.size(); ++i)
|
||||
{
|
||||
if (min_y <= ys[i] && ys[i] <= max_y)
|
||||
{
|
||||
PreciseReal y_cord_dataspace = ys[i];
|
||||
PreciseReal y = y_scale(y_cord_dataspace);
|
||||
fs << "<line x1='0' y1='" << y << "' x2='" << graph_width
|
||||
<< "' y2='" << y
|
||||
<< "' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />\n";
|
||||
|
||||
fs << "<text x='" << -margin_left/2 << "' y='" << y - 3
|
||||
<< "' font-family='times' font-size='10' fill='" << plot.font_color_ << "' transform='rotate(-90 "
|
||||
<< -margin_left/2 + 7 << " " << y << ")'>"
|
||||
<< std::setprecision(4) << y_cord_dataspace << "</text>\n";
|
||||
}
|
||||
}
|
||||
for (int i = 1; i <= plot.vertical_lines_; ++i)
|
||||
{
|
||||
CoarseReal x_cord_dataspace = plot.a_ + ((plot.b_ - plot.a_)*i)/plot.vertical_lines_;
|
||||
CoarseReal x = x_scale(x_cord_dataspace);
|
||||
fs << "<line x1='" << x << "' y1='0' x2='" << x
|
||||
<< "' y2='" << graph_height
|
||||
<< "' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />\n";
|
||||
|
||||
fs << "<text x='" << x - 10 << "' y='" << graph_height + 10
|
||||
<< "' font-family='times' font-size='10' fill='" << plot.font_color_ << "'>"
|
||||
<< std::setprecision(4) << x_cord_dataspace << "</text>\n";
|
||||
}
|
||||
}
|
||||
|
||||
int color_idx = 0;
|
||||
for (auto const & ulp : plot.ulp_list_)
|
||||
{
|
||||
std::string color = plot.colors_[color_idx++];
|
||||
for (size_t j = 0; j < ulp.size(); ++j)
|
||||
{
|
||||
if (isnan(ulp[j]))
|
||||
{
|
||||
continue;
|
||||
}
|
||||
if (plot.clip_ > 0 && abs(ulp[j]) > plot.clip_)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
CoarseReal x = x_scale(plot.coarse_abscissas_[j]);
|
||||
PreciseReal y = y_scale(ulp[j]);
|
||||
fs << "<circle cx='" << x << "' cy='" << y << "' r='1' fill='" << color << "'/>\n";
|
||||
}
|
||||
}
|
||||
|
||||
if (plot.ulp_envelope_)
|
||||
{
|
||||
std::string close_path = "' stroke='" + plot.envelope_color_ + "' stroke-width='1' fill='none'></path>\n";
|
||||
size_t jstart = 0;
|
||||
while (plot.cond_[jstart] > max_y)
|
||||
{
|
||||
++jstart;
|
||||
if (jstart >= plot.cond_.size())
|
||||
{
|
||||
goto done;
|
||||
}
|
||||
}
|
||||
|
||||
size_t jmin = jstart;
|
||||
new_top_path:
|
||||
if (jmin >= plot.cond_.size())
|
||||
{
|
||||
goto start_bottom_paths;
|
||||
}
|
||||
fs << "<path d='M" << x_scale(plot.coarse_abscissas_[jmin]) << " " << y_scale(plot.cond_[jmin]);
|
||||
|
||||
for (size_t j = jmin + 1; j < plot.coarse_abscissas_.size(); ++j)
|
||||
{
|
||||
bool bad = isnan(plot.cond_[j]) || (plot.cond_[j] > max_y);
|
||||
if (bad)
|
||||
{
|
||||
++j;
|
||||
while ( (j < plot.coarse_abscissas_.size() - 2) && bad)
|
||||
{
|
||||
bad = isnan(plot.cond_[j]) || (plot.cond_[j] > max_y);
|
||||
++j;
|
||||
}
|
||||
jmin = j;
|
||||
fs << close_path;
|
||||
goto new_top_path;
|
||||
}
|
||||
|
||||
CoarseReal t = x_scale(plot.coarse_abscissas_[j]);
|
||||
PreciseReal y = y_scale(plot.cond_[j]);
|
||||
fs << " L" << t << " " << y;
|
||||
}
|
||||
fs << close_path;
|
||||
start_bottom_paths:
|
||||
jmin = jstart;
|
||||
new_bottom_path:
|
||||
if (jmin >= plot.cond_.size())
|
||||
{
|
||||
goto done;
|
||||
}
|
||||
fs << "<path d='M" << x_scale(plot.coarse_abscissas_[jmin]) << " " << y_scale(-plot.cond_[jmin]);
|
||||
|
||||
for (size_t j = jmin + 1; j < plot.coarse_abscissas_.size(); ++j)
|
||||
{
|
||||
bool bad = isnan(plot.cond_[j]) || (-plot.cond_[j] < min_y);
|
||||
if (bad)
|
||||
{
|
||||
++j;
|
||||
while ( (j < plot.coarse_abscissas_.size() - 2) && bad)
|
||||
{
|
||||
bad = isnan(plot.cond_[j]) || (-plot.cond_[j] < min_y);
|
||||
++j;
|
||||
}
|
||||
jmin = j;
|
||||
fs << close_path;
|
||||
goto new_bottom_path;
|
||||
}
|
||||
CoarseReal t = x_scale(plot.coarse_abscissas_[j]);
|
||||
PreciseReal y = y_scale(-plot.cond_[j]);
|
||||
fs << " L" << t << " " << y;
|
||||
}
|
||||
fs << close_path;
|
||||
}
|
||||
done:
|
||||
fs << "</g>\n"
|
||||
<< "</svg>\n";
|
||||
return fs;
|
||||
}
|
||||
|
||||
private:
|
||||
std::vector<PreciseReal> precise_abscissas_;
|
||||
std::vector<CoarseReal> coarse_abscissas_;
|
||||
std::vector<PreciseReal> precise_ordinates_;
|
||||
std::vector<PreciseReal> cond_;
|
||||
std::list<std::vector<CoarseReal>> ulp_list_;
|
||||
std::vector<std::string> colors_;
|
||||
CoarseReal a_;
|
||||
CoarseReal b_;
|
||||
PreciseReal clip_;
|
||||
int width_;
|
||||
std::string envelope_color_;
|
||||
bool ulp_envelope_;
|
||||
int horizontal_lines_;
|
||||
int vertical_lines_;
|
||||
std::string title_;
|
||||
std::string background_color_;
|
||||
std::string font_color_;
|
||||
};
|
||||
|
||||
template<class F, typename PreciseReal, typename CoarseReal>
|
||||
ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::envelope_color(std::string const & color)
|
||||
{
|
||||
envelope_color_ = color;
|
||||
return *this;
|
||||
}
|
||||
|
||||
template<class F, typename PreciseReal, typename CoarseReal>
|
||||
ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::clip(PreciseReal clip)
|
||||
{
|
||||
clip_ = clip;
|
||||
return *this;
|
||||
}
|
||||
|
||||
template<class F, typename PreciseReal, typename CoarseReal>
|
||||
ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::width(int width)
|
||||
{
|
||||
width_ = width;
|
||||
return *this;
|
||||
}
|
||||
|
||||
template<class F, typename PreciseReal, typename CoarseReal>
|
||||
ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::horizontal_lines(int horizontal_lines)
|
||||
{
|
||||
horizontal_lines_ = horizontal_lines;
|
||||
return *this;
|
||||
}
|
||||
|
||||
template<class F, typename PreciseReal, typename CoarseReal>
|
||||
ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::vertical_lines(int vertical_lines)
|
||||
{
|
||||
vertical_lines_ = vertical_lines;
|
||||
return *this;
|
||||
}
|
||||
|
||||
template<class F, typename PreciseReal, typename CoarseReal>
|
||||
ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::title(std::string const & title)
|
||||
{
|
||||
title_ = title;
|
||||
return *this;
|
||||
}
|
||||
|
||||
template<class F, typename PreciseReal, typename CoarseReal>
|
||||
ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::background_color(std::string const & background_color)
|
||||
{
|
||||
background_color_ = background_color;
|
||||
return *this;
|
||||
}
|
||||
|
||||
template<class F, typename PreciseReal, typename CoarseReal>
|
||||
ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::font_color(std::string const & font_color)
|
||||
{
|
||||
font_color_ = font_color;
|
||||
return *this;
|
||||
}
|
||||
|
||||
template<class F, typename PreciseReal, typename CoarseReal>
|
||||
ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::ulp_envelope(bool write_ulp_envelope)
|
||||
{
|
||||
ulp_envelope_ = write_ulp_envelope;
|
||||
return *this;
|
||||
}
|
||||
|
||||
template<class F, typename PreciseReal, typename CoarseReal>
|
||||
void ulps_plot<F, PreciseReal, CoarseReal>::write(std::string const & filename) const
|
||||
{
|
||||
if (!boost::algorithm::ends_with(filename, ".svg"))
|
||||
{
|
||||
throw std::logic_error("Only svg files are supported at this time.");
|
||||
}
|
||||
std::ofstream fs(filename);
|
||||
fs << *this;
|
||||
fs.close();
|
||||
}
|
||||
|
||||
|
||||
template<class F, typename PreciseReal, typename CoarseReal>
|
||||
ulps_plot<F, PreciseReal, CoarseReal>::ulps_plot(F hi_acc_impl, CoarseReal a, CoarseReal b,
|
||||
size_t samples, bool perturb_abscissas, int random_seed)
|
||||
{
|
||||
static_assert(sizeof(PreciseReal) >= sizeof(CoarseReal), "PreciseReal must have larger size than CoarseReal");
|
||||
if (samples < 10)
|
||||
{
|
||||
throw std::domain_error("Must have at least 10 samples, samples = " + std::to_string(samples));
|
||||
}
|
||||
if (b <= a)
|
||||
{
|
||||
throw std::domain_error("On interval [a,b], b > a is required.");
|
||||
}
|
||||
a_ = a;
|
||||
b_ = b;
|
||||
|
||||
std::mt19937_64 gen;
|
||||
if (random_seed == -1)
|
||||
{
|
||||
std::random_device rd;
|
||||
gen.seed(rd());
|
||||
}
|
||||
// Boost's uniform_real_distribution can generate quad and multiprecision random numbers; std's cannot:
|
||||
boost::random::uniform_real_distribution<PreciseReal> dis(a, b);
|
||||
precise_abscissas_.resize(samples);
|
||||
coarse_abscissas_.resize(samples);
|
||||
|
||||
if (perturb_abscissas)
|
||||
{
|
||||
for(size_t i = 0; i < samples; ++i)
|
||||
{
|
||||
precise_abscissas_[i] = dis(gen);
|
||||
}
|
||||
std::sort(precise_abscissas_.begin(), precise_abscissas_.end());
|
||||
for (size_t i = 0; i < samples; ++i)
|
||||
{
|
||||
coarse_abscissas_[i] = static_cast<CoarseReal>(precise_abscissas_[i]);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
for(size_t i = 0; i < samples; ++i)
|
||||
{
|
||||
coarse_abscissas_[i] = static_cast<CoarseReal>(dis(gen));
|
||||
}
|
||||
std::sort(coarse_abscissas_.begin(), coarse_abscissas_.end());
|
||||
for (size_t i = 0; i < samples; ++i)
|
||||
{
|
||||
precise_abscissas_[i] = coarse_abscissas_[i];
|
||||
}
|
||||
}
|
||||
|
||||
precise_ordinates_.resize(samples);
|
||||
for (size_t i = 0; i < samples; ++i)
|
||||
{
|
||||
precise_ordinates_[i] = hi_acc_impl(precise_abscissas_[i]);
|
||||
}
|
||||
|
||||
cond_.resize(samples, std::numeric_limits<PreciseReal>::quiet_NaN());
|
||||
for (size_t i = 0 ; i < samples; ++i)
|
||||
{
|
||||
PreciseReal y = precise_ordinates_[i];
|
||||
if (y != 0)
|
||||
{
|
||||
// Maybe cond_ is badly names; should it be half_cond_?
|
||||
cond_[i] = boost::math::tools::evaluation_condition_number(hi_acc_impl, precise_abscissas_[i])/2;
|
||||
// Half-ULP accuracy is the correctly rounded result, so make sure the envelop doesn't go below this:
|
||||
if (cond_[i] < 0.5)
|
||||
{
|
||||
cond_[i] = 0.5;
|
||||
}
|
||||
}
|
||||
// else leave it as nan.
|
||||
}
|
||||
clip_ = -1;
|
||||
width_ = 1100;
|
||||
envelope_color_ = "chartreuse";
|
||||
ulp_envelope_ = true;
|
||||
horizontal_lines_ = 8;
|
||||
vertical_lines_ = 10;
|
||||
title_ = "";
|
||||
background_color_ = "black";
|
||||
font_color_ = "white";
|
||||
}
|
||||
|
||||
template<class F, typename PreciseReal, typename CoarseReal>
|
||||
template<class G>
|
||||
ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::add_fn(G g, std::string const & color)
|
||||
{
|
||||
using std::abs;
|
||||
size_t samples = precise_abscissas_.size();
|
||||
std::vector<CoarseReal> ulps(samples);
|
||||
for (size_t i = 0; i < samples; ++i)
|
||||
{
|
||||
PreciseReal y_hi_acc = precise_ordinates_[i];
|
||||
PreciseReal y_lo_acc = g(coarse_abscissas_[i]);
|
||||
PreciseReal absy = abs(y_hi_acc);
|
||||
PreciseReal dist = nextafter(static_cast<CoarseReal>(absy), std::numeric_limits<CoarseReal>::max()) - static_cast<CoarseReal>(absy);
|
||||
ulps[i] = static_cast<CoarseReal>((y_lo_acc - y_hi_acc)/dist);
|
||||
}
|
||||
ulp_list_.emplace_back(ulps);
|
||||
colors_.emplace_back(color);
|
||||
return *this;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
} // namespace boost::math::tools
|
||||
#endif
|
||||
57
reporting/performance/chebyshev_clenshaw.cpp
Normal file
57
reporting/performance/chebyshev_clenshaw.cpp
Normal file
@@ -0,0 +1,57 @@
|
||||
#include <vector>
|
||||
#include <random>
|
||||
#include <benchmark/benchmark.h>
|
||||
#include <boost/math/special_functions/chebyshev.hpp>
|
||||
|
||||
|
||||
template<class Real>
|
||||
void ChebyshevClenshaw(benchmark::State& state)
|
||||
{
|
||||
std::vector<Real> v(state.range(0));
|
||||
std::random_device rd;
|
||||
std::mt19937_64 mt(rd());
|
||||
std::uniform_real_distribution<Real> unif(-1,1);
|
||||
for (size_t i = 0; i < v.size(); ++i)
|
||||
{
|
||||
v[i] = unif(mt);
|
||||
}
|
||||
|
||||
using boost::math::chebyshev_clenshaw_recurrence;
|
||||
Real x = unif(mt);
|
||||
for (auto _ : state)
|
||||
{
|
||||
benchmark::DoNotOptimize(chebyshev_clenshaw_recurrence(v.data(), v.size(), x));
|
||||
}
|
||||
state.SetComplexityN(state.range(0));
|
||||
}
|
||||
|
||||
template<class Real>
|
||||
void TranslatedChebyshevClenshaw(benchmark::State& state)
|
||||
{
|
||||
std::vector<Real> v(state.range(0));
|
||||
std::random_device rd;
|
||||
std::mt19937_64 mt(rd());
|
||||
std::uniform_real_distribution<Real> unif(-1,1);
|
||||
for (size_t i = 0; i < v.size(); ++i)
|
||||
{
|
||||
v[i] = unif(mt);
|
||||
}
|
||||
|
||||
using boost::math::detail::unchecked_chebyshev_clenshaw_recurrence;
|
||||
Real x = unif(mt);
|
||||
Real a = -2;
|
||||
Real b = 5;
|
||||
for (auto _ : state)
|
||||
{
|
||||
benchmark::DoNotOptimize(unchecked_chebyshev_clenshaw_recurrence(v.data(), v.size(), a, b, x));
|
||||
}
|
||||
state.SetComplexityN(state.range(0));
|
||||
}
|
||||
|
||||
|
||||
BENCHMARK_TEMPLATE(TranslatedChebyshevClenshaw, double)->RangeMultiplier(2)->Range(1<<1, 1<<22)->Complexity(benchmark::oN);
|
||||
BENCHMARK_TEMPLATE(ChebyshevClenshaw, double)->RangeMultiplier(2)->Range(1<<1, 1<<22)->Complexity(benchmark::oN);
|
||||
|
||||
|
||||
|
||||
BENCHMARK_MAIN();
|
||||
@@ -464,11 +464,11 @@ test-suite special_fun :
|
||||
[ run test_lambert_w_derivative.cpp ../../test/build//boost_unit_test_framework : : : <define>BOOST_MATH_TEST_MULTIPRECISION [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <define>BOOST_MATH_TEST_FLOAT128 <linkflags>"-Bstatic -lquadmath -Bdynamic" ] ]
|
||||
|
||||
[ run test_legendre.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework : : : [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>"-Bstatic -lquadmath -Bdynamic" ] ]
|
||||
[ run chebyshev_test.cpp : : : [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>"-Bstatic -lquadmath -Bdynamic" ] ]
|
||||
[ run chebyshev_transform_test.cpp ../config//fftw3f : : : <define>TEST1 [ requires cxx11_auto_declarations cxx11_range_based_for ] [ check-target-builds ../config//has_fftw3 "libfftw3" : : <build>no ] : chebyshev_transform_test_1 ]
|
||||
[ run chebyshev_transform_test.cpp ../config//fftw3 : : : <define>TEST2 [ requires cxx11_auto_declarations cxx11_range_based_for ] [ check-target-builds ../config//has_fftw3 "libfftw3" : : <build>no ] : chebyshev_transform_test_2 ]
|
||||
[ run chebyshev_transform_test.cpp ../config//fftw3l : : : <define>TEST3 [ requires cxx11_auto_declarations cxx11_range_based_for ] [ check-target-builds ../config//has_fftw3 "libfftw3" : : <build>no ] : chebyshev_transform_test_3 ]
|
||||
[ run chebyshev_transform_test.cpp ../config//fftw3q ../config//quadmath : : : <define>TEST4 [ requires cxx11_auto_declarations cxx11_range_based_for ] [ check-target-builds ../config//has_fftw3 "libfftw3" : : <build>no ] [ check-target-builds ../config//has_float128 "__float128" : : <build>no ] : chebyshev_transform_test_4 ]
|
||||
[ run chebyshev_test.cpp : : : [ requires cxx11_inline_namespaces cxx11_unified_initialization_syntax cxx11_hdr_tuple cxx11_smart_ptr cxx11_defaulted_functions cxx11_auto_declarations cxx11_range_based_for cxx11_constexpr ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>"-Bstatic -lquadmath -Bdynamic" ] ]
|
||||
[ run chebyshev_transform_test.cpp ../config//fftw3f : : : <define>TEST1 [ requires cxx11_smart_ptr cxx11_defaulted_functions cxx11_auto_declarations cxx11_range_based_for ] [ check-target-builds ../config//has_fftw3 "libfftw3" : : <build>no ] : chebyshev_transform_test_1 ]
|
||||
[ run chebyshev_transform_test.cpp ../config//fftw3 : : : <define>TEST2 [ requires cxx11_smart_ptr cxx11_defaulted_functions cxx11_auto_declarations cxx11_range_based_for ] [ check-target-builds ../config//has_fftw3 "libfftw3" : : <build>no ] : chebyshev_transform_test_2 ]
|
||||
[ run chebyshev_transform_test.cpp ../config//fftw3l : : : <define>TEST3 [ requires cxx11_smart_ptr cxx11_defaulted_functions cxx11_auto_declarations cxx11_range_based_for ] [ check-target-builds ../config//has_fftw3 "libfftw3" : : <build>no ] : chebyshev_transform_test_3 ]
|
||||
[ run chebyshev_transform_test.cpp ../config//fftw3q ../config//quadmath : : : <define>TEST4 [ requires cxx11_smart_ptr cxx11_defaulted_functions cxx11_auto_declarations cxx11_range_based_for ] [ check-target-builds ../config//has_fftw3 "libfftw3" : : <build>no ] [ check-target-builds ../config//has_float128 "__float128" : : <build>no ] : chebyshev_transform_test_4 ]
|
||||
|
||||
[ run cardinal_trigonometric_test.cpp ../config//fftw3f : : : <define>TEST1 [ requires cxx11_auto_declarations cxx11_range_based_for ] [ check-target-builds ../config//has_fftw3 "libfftw3" : : <build>no ] : cardinal_trigonometric_test_1 ]
|
||||
[ run cardinal_trigonometric_test.cpp ../config//fftw3 : : : <define>TEST2 [ requires cxx11_auto_declarations cxx11_range_based_for ] [ check-target-builds ../config//has_fftw3 "libfftw3" : : <build>no ] : cardinal_trigonometric_test_2 ]
|
||||
|
||||
@@ -5,20 +5,16 @@
|
||||
* LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||||
*/
|
||||
|
||||
#define BOOST_TEST_MODULE chebyshev_test
|
||||
|
||||
#include <random>
|
||||
#include <iostream>
|
||||
#include <boost/type_index.hpp>
|
||||
#include <boost/test/included/unit_test.hpp>
|
||||
#include <boost/test/tools/floating_point_comparison.hpp>
|
||||
#include <boost/math/special_functions/chebyshev.hpp>
|
||||
#include <boost/math/special_functions/sinc.hpp>
|
||||
#include <boost/multiprecision/cpp_bin_float.hpp>
|
||||
#include <boost/multiprecision/cpp_dec_float.hpp>
|
||||
#include <boost/array.hpp>
|
||||
#include "math_unit_test.hpp"
|
||||
|
||||
using boost::multiprecision::cpp_bin_float_quad;
|
||||
using boost::multiprecision::cpp_bin_float_50;
|
||||
using boost::multiprecision::cpp_bin_float_100;
|
||||
using boost::math::chebyshev_t;
|
||||
using boost::math::chebyshev_t_prime;
|
||||
using boost::math::chebyshev_u;
|
||||
@@ -29,28 +25,25 @@ void test_polynomials()
|
||||
std::cout << "Testing explicit polynomial representations of the Chebyshev polynomials on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
|
||||
|
||||
Real x = -2;
|
||||
Real tol = 400*std::numeric_limits<Real>::epsilon();
|
||||
if (tol > std::numeric_limits<float>::epsilon())
|
||||
tol *= 10; // float results have much larger error rates.
|
||||
Real tol = 32*std::numeric_limits<Real>::epsilon();
|
||||
while (x < 2)
|
||||
{
|
||||
BOOST_CHECK_CLOSE_FRACTION(chebyshev_t(0, x), Real(1), tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(chebyshev_t(1, x), x, tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(chebyshev_t(2, x), 2*x*x - 1, tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(chebyshev_t(3, x), x*(4*x*x-3), tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(chebyshev_t(4, x), 8*x*x*(x*x - 1) + 1, tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(chebyshev_t(5, x), x*(16*x*x*x*x - 20*x*x + 5), tol);
|
||||
CHECK_MOLLIFIED_CLOSE(chebyshev_t(0, x), Real(1), tol);
|
||||
CHECK_MOLLIFIED_CLOSE(chebyshev_t(1, x), x, tol);
|
||||
CHECK_MOLLIFIED_CLOSE(chebyshev_t(2, x), 2*x*x - 1, tol);
|
||||
CHECK_MOLLIFIED_CLOSE(chebyshev_t(3, x), x*(4*x*x-3), tol);
|
||||
CHECK_MOLLIFIED_CLOSE(chebyshev_t(4, x), 8*x*x*(x*x - 1) + 1, tol);
|
||||
CHECK_MOLLIFIED_CLOSE(chebyshev_t(5, x), x*(16*x*x*x*x - 20*x*x + 5), tol);
|
||||
x += 1/static_cast<Real>(1<<7);
|
||||
}
|
||||
|
||||
x = -2;
|
||||
tol = 10*tol;
|
||||
while (x < 2)
|
||||
{
|
||||
BOOST_CHECK_CLOSE_FRACTION(chebyshev_u(0, x), Real(1), tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(chebyshev_u(1, x), 2*x, tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(chebyshev_u(2, x), 4*x*x - 1, tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(chebyshev_u(3, x), 4*x*(2*x*x - 1), tol);
|
||||
CHECK_MOLLIFIED_CLOSE(chebyshev_u(0, x), Real(1), tol);
|
||||
CHECK_MOLLIFIED_CLOSE(chebyshev_u(1, x), 2*x, tol);
|
||||
CHECK_MOLLIFIED_CLOSE(chebyshev_u(2, x), 4*x*x - 1, tol);
|
||||
CHECK_MOLLIFIED_CLOSE(chebyshev_u(3, x), 4*x*(2*x*x - 1), tol);
|
||||
x += 1/static_cast<Real>(1<<7);
|
||||
}
|
||||
}
|
||||
@@ -62,14 +55,14 @@ void test_derivatives()
|
||||
std::cout << "Testing explicit polynomial representations of the Chebyshev polynomial derivatives on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
|
||||
|
||||
Real x = -2;
|
||||
Real tol = 1000*std::numeric_limits<Real>::epsilon();
|
||||
Real tol = 4*std::numeric_limits<Real>::epsilon();
|
||||
while (x < 2)
|
||||
{
|
||||
BOOST_CHECK_CLOSE_FRACTION(chebyshev_t_prime(0, x), Real(0), tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(chebyshev_t_prime(1, x), Real(1), tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(chebyshev_t_prime(2, x), 4*x, tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(chebyshev_t_prime(3, x), 3*(4*x*x - 1), tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(chebyshev_t_prime(4, x), 16*x*(2*x*x - 1), tol);
|
||||
CHECK_MOLLIFIED_CLOSE(chebyshev_t_prime(0, x), Real(0), tol);
|
||||
CHECK_MOLLIFIED_CLOSE(chebyshev_t_prime(1, x), Real(1), tol);
|
||||
CHECK_MOLLIFIED_CLOSE(chebyshev_t_prime(2, x), 4*x, tol);
|
||||
CHECK_MOLLIFIED_CLOSE(chebyshev_t_prime(3, x), 3*(4*x*x - 1), tol);
|
||||
CHECK_MOLLIFIED_CLOSE(chebyshev_t_prime(4, x), 16*x*(2*x*x - 1), tol);
|
||||
// This one makes the tolerance have to grow too large; the Chebyshev recurrence is more stable than naive polynomial evaluation anyway.
|
||||
//BOOST_CHECK_CLOSE_FRACTION(chebyshev_t_prime(5, x), 5*(4*x*x*(4*x*x - 3) + 1), tol);
|
||||
x += 1/static_cast<Real>(1<<7);
|
||||
@@ -93,48 +86,93 @@ void test_clenshaw_recurrence()
|
||||
boost::array<Real, 7> c6 = { {0, 0, 0, 0, 0, 0, 1} };
|
||||
|
||||
Real x = -1;
|
||||
Real tol = 10*std::numeric_limits<Real>::epsilon();
|
||||
if (tol > std::numeric_limits<float>::epsilon())
|
||||
tol *= 100; // float results have much larger error rates.
|
||||
// It's not clear from this test which one is more accurate; higher precision cast testing is required, and is done elsewhere:
|
||||
int ulps = 50;
|
||||
while (x <= 1)
|
||||
{
|
||||
Real y = chebyshev_clenshaw_recurrence(c0.data(), c0.size(), x);
|
||||
BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(0, x), tol);
|
||||
Real expected = chebyshev_t(0, x);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
y = chebyshev_clenshaw_recurrence(c0.data(), c0.size(), Real(-1), Real(1), x);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
|
||||
y = chebyshev_clenshaw_recurrence(c01.data(), c01.size(), x);
|
||||
BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(0, x), tol);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
y = chebyshev_clenshaw_recurrence(c01.data(), c01.size(), Real(-1), Real(1), x);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
|
||||
y = chebyshev_clenshaw_recurrence(c02.data(), c02.size(), x);
|
||||
BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(0, x), tol);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
y = chebyshev_clenshaw_recurrence(c02.data(), c02.size(), Real(-1), Real(1), x);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
|
||||
expected = chebyshev_t(1,x);
|
||||
y = chebyshev_clenshaw_recurrence(c1.data(), c1.size(), x);
|
||||
BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(1, x), tol);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
y = chebyshev_clenshaw_recurrence(c1.data(), c1.size(), Real(-1), Real(1), x);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
|
||||
expected = chebyshev_t(2, x);
|
||||
y = chebyshev_clenshaw_recurrence(c2.data(), c2.size(), x);
|
||||
BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(2, x), tol);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
y = chebyshev_clenshaw_recurrence(c2.data(), c2.size(), Real(-1), Real(1), x);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
|
||||
expected = chebyshev_t(3, x);
|
||||
y = chebyshev_clenshaw_recurrence(c3.data(), c3.size(), x);
|
||||
BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(3, x), tol);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
y = chebyshev_clenshaw_recurrence(c3.data(), c3.size(), Real(-1), Real(1), x);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
|
||||
expected = chebyshev_t(4, x);
|
||||
y = chebyshev_clenshaw_recurrence(c4.data(), c4.size(), x);
|
||||
BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(4, x), tol);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
y = chebyshev_clenshaw_recurrence(c4.data(), c4.size(), Real(-1), Real(1), x);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
|
||||
expected = chebyshev_t(5, x);
|
||||
y = chebyshev_clenshaw_recurrence(c5.data(), c5.size(), x);
|
||||
BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(5, x), tol);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
y = chebyshev_clenshaw_recurrence(c5.data(), c5.size(), Real(-1), Real(1), x);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
|
||||
expected = chebyshev_t(6, x);
|
||||
y = chebyshev_clenshaw_recurrence(c6.data(), c6.size(), x);
|
||||
BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(6, x), tol);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
y = chebyshev_clenshaw_recurrence(c6.data(), c6.size(), Real(-1), Real(1), x);
|
||||
CHECK_ULP_CLOSE(expected, y, ulps);
|
||||
|
||||
x += static_cast<Real>(1)/static_cast<Real>(1 << 7);
|
||||
}
|
||||
}
|
||||
|
||||
BOOST_AUTO_TEST_CASE(chebyshev_test)
|
||||
template<typename Real>
|
||||
void test_translated_clenshaw_recurrence()
|
||||
{
|
||||
test_clenshaw_recurrence<float>();
|
||||
test_clenshaw_recurrence<double>();
|
||||
test_clenshaw_recurrence<long double>();
|
||||
using boost::math::chebyshev_clenshaw_recurrence;
|
||||
std::mt19937_64 mt(123242);
|
||||
std::uniform_real_distribution<Real> dis(-1,1);
|
||||
|
||||
std::vector<Real> c(32);
|
||||
for (auto & d : c) {
|
||||
d = dis(mt);
|
||||
}
|
||||
int samples = 0;
|
||||
while (samples++ < 250) {
|
||||
Real x = dis(mt);
|
||||
Real expected = chebyshev_clenshaw_recurrence(c.data(), c.size(), x);
|
||||
// The computed value, in this case, should actually be better, since it has more information to use.
|
||||
// In any case, this test doesn't show Reinch's modification to the Clenshaw recurrence is better;
|
||||
// it shows they are doing roughly the same thing.
|
||||
Real computed = chebyshev_clenshaw_recurrence(c.data(), c.size(), Real(-1), Real(1), x);
|
||||
if (!CHECK_ULP_CLOSE(expected, computed, 1000)) {
|
||||
std::cerr << " Problem occured at x = " << x << "\n";
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
test_polynomials<float>();
|
||||
test_polynomials<double>();
|
||||
test_polynomials<long double>();
|
||||
@@ -144,4 +182,11 @@ BOOST_AUTO_TEST_CASE(chebyshev_test)
|
||||
test_derivatives<double>();
|
||||
test_derivatives<long double>();
|
||||
test_derivatives<cpp_bin_float_quad>();
|
||||
|
||||
test_clenshaw_recurrence<float>();
|
||||
test_clenshaw_recurrence<double>();
|
||||
test_clenshaw_recurrence<long double>();
|
||||
|
||||
test_translated_clenshaw_recurrence<double>();
|
||||
return boost::math::test::report_errors();
|
||||
}
|
||||
|
||||
@@ -4,17 +4,11 @@
|
||||
* Boost Software License, Version 1.0. (See accompanying file
|
||||
* LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||||
*/
|
||||
#define BOOST_TEST_MODULE chebyshev_transform_test
|
||||
|
||||
#include <boost/cstdfloat.hpp>
|
||||
#include "math_unit_test.hpp"
|
||||
#include <boost/type_index.hpp>
|
||||
#include <boost/test/included/unit_test.hpp>
|
||||
#include <boost/test/tools/floating_point_comparison.hpp>
|
||||
#include <boost/math/special_functions/chebyshev.hpp>
|
||||
#include <boost/math/special_functions/chebyshev_transform.hpp>
|
||||
#include <boost/math/special_functions/sinc.hpp>
|
||||
#include <boost/multiprecision/cpp_bin_float.hpp>
|
||||
#include <boost/multiprecision/cpp_dec_float.hpp>
|
||||
|
||||
#if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4)
|
||||
# define TEST1
|
||||
@@ -23,9 +17,6 @@
|
||||
# define TEST4
|
||||
#endif
|
||||
|
||||
using boost::multiprecision::cpp_bin_float_quad;
|
||||
using boost::multiprecision::cpp_bin_float_50;
|
||||
using boost::multiprecision::cpp_bin_float_100;
|
||||
using boost::math::chebyshev_t;
|
||||
using boost::math::chebyshev_t_prime;
|
||||
using boost::math::chebyshev_u;
|
||||
@@ -41,8 +32,8 @@ void test_sin_chebyshev_transform()
|
||||
using std::cos;
|
||||
using std::abs;
|
||||
|
||||
Real tol = 10*std::numeric_limits<Real>::epsilon();
|
||||
auto f = [](Real x) { return sin(x); };
|
||||
Real tol = std::numeric_limits<Real>::epsilon();
|
||||
auto f = [](Real x)->Real { return sin(x); };
|
||||
Real a = 0;
|
||||
Real b = 1;
|
||||
chebyshev_transform<Real> cheb(f, a, b, tol);
|
||||
@@ -52,29 +43,14 @@ void test_sin_chebyshev_transform()
|
||||
{
|
||||
Real s = sin(x);
|
||||
Real c = cos(x);
|
||||
if (abs(s) < tol)
|
||||
{
|
||||
BOOST_CHECK_SMALL(cheb(x), 100*tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(c, cheb.prime(x), 100*tol);
|
||||
}
|
||||
else
|
||||
{
|
||||
BOOST_CHECK_CLOSE_FRACTION(s, cheb(x), 100*tol);
|
||||
if (abs(c) < tol)
|
||||
{
|
||||
BOOST_CHECK_SMALL(cheb.prime(x), 100*tol);
|
||||
}
|
||||
else
|
||||
{
|
||||
BOOST_CHECK_CLOSE_FRACTION(c, cheb.prime(x), 100*tol);
|
||||
}
|
||||
}
|
||||
CHECK_ABSOLUTE_ERROR(s, cheb(x), tol);
|
||||
CHECK_ABSOLUTE_ERROR(c, cheb.prime(x), 150*tol);
|
||||
x += static_cast<Real>(1)/static_cast<Real>(1 << 7);
|
||||
}
|
||||
|
||||
Real Q = cheb.integrate();
|
||||
|
||||
BOOST_CHECK_CLOSE_FRACTION(1 - cos(static_cast<Real>(1)), Q, 100*tol);
|
||||
CHECK_ABSOLUTE_ERROR(1 - cos(static_cast<Real>(1)), Q, 100*tol);
|
||||
}
|
||||
|
||||
|
||||
@@ -88,7 +64,7 @@ void test_sinc_chebyshev_transform()
|
||||
using boost::math::chebyshev_transform;
|
||||
using boost::math::constants::half_pi;
|
||||
|
||||
Real tol = 500*std::numeric_limits<Real>::epsilon();
|
||||
Real tol = 100*std::numeric_limits<Real>::epsilon();
|
||||
auto f = [](Real x) { return boost::math::sinc_pi(x); };
|
||||
Real a = 0;
|
||||
Real b = 1;
|
||||
@@ -100,30 +76,16 @@ void test_sinc_chebyshev_transform()
|
||||
Real s = sinc_pi(x);
|
||||
Real ds = (cos(x)-sinc_pi(x))/x;
|
||||
if (x == 0) { ds = 0; }
|
||||
if (s < tol)
|
||||
{
|
||||
BOOST_CHECK_SMALL(cheb(x), tol);
|
||||
}
|
||||
else
|
||||
{
|
||||
BOOST_CHECK_CLOSE_FRACTION(s, cheb(x), tol);
|
||||
}
|
||||
|
||||
if (abs(ds) < tol)
|
||||
{
|
||||
BOOST_CHECK_SMALL(cheb.prime(x), 5 * tol);
|
||||
}
|
||||
else
|
||||
{
|
||||
BOOST_CHECK_CLOSE_FRACTION(ds, cheb.prime(x), 300*tol);
|
||||
}
|
||||
CHECK_ABSOLUTE_ERROR(s, cheb(x), tol);
|
||||
CHECK_ABSOLUTE_ERROR(ds, cheb.prime(x), 10*tol);
|
||||
x += static_cast<Real>(1)/static_cast<Real>(1 << 7);
|
||||
}
|
||||
|
||||
Real Q = cheb.integrate();
|
||||
//NIntegrate[Sinc[x], {x, 0, 1}, WorkingPrecision -> 200, AccuracyGoal -> 150, PrecisionGoal -> 150, MaxRecursion -> 150]
|
||||
Real Q_exp = boost::lexical_cast<Real>("0.94608307036718301494135331382317965781233795473811179047145477356668");
|
||||
BOOST_CHECK_CLOSE_FRACTION(Q_exp, Q, tol);
|
||||
CHECK_ABSOLUTE_ERROR(Q_exp, Q, tol);
|
||||
}
|
||||
|
||||
|
||||
@@ -133,6 +95,8 @@ template<class Real>
|
||||
void test_atap_examples()
|
||||
{
|
||||
using std::sin;
|
||||
using std::exp;
|
||||
using std::sqrt;
|
||||
using boost::math::constants::half;
|
||||
using boost::math::sinc_pi;
|
||||
using boost::math::chebyshev_transform;
|
||||
@@ -144,33 +108,37 @@ void test_atap_examples()
|
||||
Real u = (0 < s) - (s < 0);
|
||||
return t + u; };
|
||||
|
||||
auto f3 = [](Real x) { return sin(6*x) + sin(60*exp(x)); };
|
||||
|
||||
auto f4 = [](Real x) { return 1/(1+1000*(x+half<Real>())*(x+half<Real>())) + 1/sqrt(1+1000*(x-.5)*(x-0.5));};
|
||||
//auto f3 = [](Real x) { return sin(6*x) + sin(60*exp(x)); };
|
||||
//auto f4 = [](Real x) { return 1/(1+1000*(x+half<Real>())*(x+half<Real>())) + 1/sqrt(1+1000*(x-Real(1)/Real(2))*(x-Real(1)/Real(2)));};
|
||||
Real a = -1;
|
||||
Real b = 1;
|
||||
chebyshev_transform<Real> cheb1(f1, a, b);
|
||||
chebyshev_transform<Real> cheb1(f1, a, b, tol);
|
||||
chebyshev_transform<Real> cheb2(f2, a, b, tol);
|
||||
//chebyshev_transform<Real> cheb3(f3, a, b, tol);
|
||||
|
||||
Real x = a;
|
||||
while (x < b)
|
||||
{
|
||||
//Real s = f1(x);
|
||||
if (sizeof(Real) == sizeof(float))
|
||||
// f1 and f2 are not differentiable; standard convergence rate theorems don't apply.
|
||||
// Basically, the max refinements are always hit; so the error is not related to the precision of the type.
|
||||
Real acceptable_error = sqrt(tol);
|
||||
Real acceptable_error_2 = 9e-4;
|
||||
if (std::is_same<Real, long double>::value)
|
||||
{
|
||||
BOOST_CHECK_CLOSE_FRACTION(f1(x), cheb1(x), 4e-3);
|
||||
acceptable_error = 1.6e-5;
|
||||
}
|
||||
else
|
||||
if (std::is_same<Real, double>::value)
|
||||
{
|
||||
BOOST_CHECK_CLOSE_FRACTION(f1(x), cheb1(x), 1.3e-5);
|
||||
acceptable_error *= 500;
|
||||
}
|
||||
BOOST_CHECK_CLOSE_FRACTION(f2(x), cheb2(x), 6e-3);
|
||||
//BOOST_CHECK_CLOSE_FRACTION(f3(x), cheb3(x), 100*tol);
|
||||
CHECK_ABSOLUTE_ERROR(f1(x), cheb1(x), acceptable_error);
|
||||
|
||||
CHECK_ABSOLUTE_ERROR(f2(x), cheb2(x), acceptable_error_2);
|
||||
x += static_cast<Real>(1)/static_cast<Real>(1 << 7);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
//Validate that the Chebyshev polynomials are well approximated by the Chebyshev transform.
|
||||
template<class Real>
|
||||
void test_chebyshev_chebyshev_transform()
|
||||
@@ -179,58 +147,44 @@ void test_chebyshev_chebyshev_transform()
|
||||
// T_0 = 1:
|
||||
auto t0 = [](Real) { return 1; };
|
||||
chebyshev_transform<Real> cheb0(t0, -1, 1);
|
||||
BOOST_CHECK_CLOSE_FRACTION(cheb0.coefficients()[0], 2, tol);
|
||||
CHECK_ABSOLUTE_ERROR(2, cheb0.coefficients()[0], tol);
|
||||
|
||||
Real x = -1;
|
||||
while (x < 1)
|
||||
{
|
||||
BOOST_CHECK_CLOSE_FRACTION(cheb0(x), 1, tol);
|
||||
BOOST_CHECK_SMALL(cheb0.prime(x), tol);
|
||||
CHECK_ABSOLUTE_ERROR(1, cheb0(x), tol);
|
||||
CHECK_ABSOLUTE_ERROR(Real(0), cheb0.prime(x), tol);
|
||||
x += static_cast<Real>(1)/static_cast<Real>(1 << 7);
|
||||
}
|
||||
|
||||
// T_1 = x:
|
||||
auto t1 = [](Real x) { return x; };
|
||||
chebyshev_transform<Real> cheb1(t1, -1, 1);
|
||||
BOOST_CHECK_CLOSE_FRACTION(cheb1.coefficients()[1], 1, tol);
|
||||
CHECK_ABSOLUTE_ERROR(Real(1), cheb1.coefficients()[1], tol);
|
||||
|
||||
x = -1;
|
||||
while (x < 1)
|
||||
{
|
||||
if (x == 0)
|
||||
{
|
||||
BOOST_CHECK_SMALL(cheb1(x), tol);
|
||||
}
|
||||
else
|
||||
{
|
||||
BOOST_CHECK_CLOSE_FRACTION(cheb1(x), x, tol);
|
||||
}
|
||||
BOOST_CHECK_CLOSE_FRACTION(cheb1.prime(x), 1, tol);
|
||||
CHECK_ABSOLUTE_ERROR(x, cheb1(x), tol);
|
||||
CHECK_ABSOLUTE_ERROR(Real(1), cheb1.prime(x), tol);
|
||||
x += static_cast<Real>(1)/static_cast<Real>(1 << 7);
|
||||
}
|
||||
|
||||
|
||||
auto t2 = [](Real x) { return 2*x*x-1; };
|
||||
chebyshev_transform<Real> cheb2(t2, -1, 1);
|
||||
BOOST_CHECK_CLOSE_FRACTION(cheb2.coefficients()[2], 1, tol);
|
||||
CHECK_ABSOLUTE_ERROR(Real(1), cheb2.coefficients()[2], tol);
|
||||
|
||||
x = -1;
|
||||
while (x < 1)
|
||||
{
|
||||
BOOST_CHECK_CLOSE_FRACTION(cheb2(x), t2(x), tol);
|
||||
if (x != 0)
|
||||
{
|
||||
BOOST_CHECK_CLOSE_FRACTION(cheb2.prime(x), 4*x, tol);
|
||||
}
|
||||
else
|
||||
{
|
||||
BOOST_CHECK_SMALL(cheb2.prime(x), tol);
|
||||
}
|
||||
CHECK_ABSOLUTE_ERROR(t2(x), cheb2(x), tol);
|
||||
CHECK_ABSOLUTE_ERROR(4*x, cheb2.prime(x), tol);
|
||||
x += static_cast<Real>(1)/static_cast<Real>(1 << 7);
|
||||
}
|
||||
}
|
||||
|
||||
BOOST_AUTO_TEST_CASE(chebyshev_transform_test)
|
||||
int main()
|
||||
{
|
||||
#ifdef TEST1
|
||||
test_chebyshev_chebyshev_transform<float>();
|
||||
@@ -258,6 +212,8 @@ BOOST_AUTO_TEST_CASE(chebyshev_transform_test)
|
||||
test_sinc_chebyshev_transform<__float128>();
|
||||
#endif
|
||||
#endif
|
||||
|
||||
return boost::math::test::report_errors();
|
||||
}
|
||||
|
||||
|
||||
|
||||
@@ -236,7 +236,7 @@ bool check_conditioned_error(Real abscissa, PreciseReal expected1, PreciseReal e
|
||||
<< " This exceeds the tolerance " << tol << "\n"
|
||||
<< std::showpos
|
||||
<< " Expected: " << std::defaultfloat << std::fixed << expected << " = " << std::scientific << expected << std::hexfloat << " = " << expected << "\n"
|
||||
<< " Computed: " << std::defaultfloat << std::fixed << computed << " = " << std::scientific << expected << std::hexfloat << " = " << computed << "\n"
|
||||
<< " Computed: " << std::defaultfloat << std::fixed << computed << " = " << std::scientific << computed << std::hexfloat << " = " << computed << "\n"
|
||||
<< " Condition number of function evaluation: " << std::noshowpos << std::defaultfloat << std::scientific << cond << " = " << std::fixed << cond << "\n"
|
||||
<< " Badness scale required to make this message go away: " << std::defaultfloat << relative_error/(cond*mu) << "\n";
|
||||
std::cerr.flags(f);
|
||||
@@ -247,6 +247,53 @@ bool check_conditioned_error(Real abscissa, PreciseReal expected1, PreciseReal e
|
||||
}
|
||||
|
||||
|
||||
template<class PreciseReal, class Real>
|
||||
bool check_absolute_error(PreciseReal expected1, Real computed, Real acceptable_error, std::string const & filename, std::string const & function, int line)
|
||||
{
|
||||
using std::max;
|
||||
using std::abs;
|
||||
using std::isnan;
|
||||
// Of course integers can be expected values, and they are exact:
|
||||
if (!std::is_integral<PreciseReal>::value) {
|
||||
BOOST_ASSERT_MSG(sizeof(PreciseReal) >= sizeof(Real),
|
||||
"The expected number must be computed in higher (or equal) precision than the number being tested.");
|
||||
BOOST_ASSERT_MSG(!isnan(expected1), "Expected value cannot be a nan.");
|
||||
}
|
||||
BOOST_ASSERT_MSG(acceptable_error > 0, "Error must be > 0.");
|
||||
|
||||
if (isnan(computed))
|
||||
{
|
||||
std::ios_base::fmtflags f(std::cerr.flags());
|
||||
std::cerr << std::setprecision(3);
|
||||
std::cerr << "\033[0;31mError at " << filename << ":" << function << ":" << line << ":\n"
|
||||
<< " \033[0m Computed value is a nan\n";
|
||||
std::cerr.flags(f);
|
||||
++detail::global_error_count;
|
||||
return false;
|
||||
}
|
||||
|
||||
Real expected = Real(expected1);
|
||||
Real error = abs(expected - computed);
|
||||
if (error > acceptable_error)
|
||||
{
|
||||
std::ios_base::fmtflags f( std::cerr.flags() );
|
||||
std::cerr << std::setprecision(3);
|
||||
std::cerr << "\033[0;31mError at " << filename << ":" << function << ":" << line << "\n";
|
||||
std::cerr << std::setprecision(std::numeric_limits<Real>::max_digits10);
|
||||
std::cerr << "\033[0m The absolute error in " << boost::core::demangle(typeid(Real).name()) << " precision is " << std::scientific << error << "\n"
|
||||
<< " This exceeds the acceptable error " << acceptable_error << "\n"
|
||||
<< std::showpos
|
||||
<< " Expected: " << std::defaultfloat << std::fixed << expected << " = " << std::scientific << expected << std::hexfloat << " = " << expected << "\n"
|
||||
<< " Computed: " << std::defaultfloat << std::fixed << computed << " = " << std::scientific << computed<< std::hexfloat << " = " << computed << "\n"
|
||||
<< " Error/Acceptable error: " << std::defaultfloat << error/acceptable_error << "\n";
|
||||
std::cerr.flags(f);
|
||||
++detail::global_error_count;
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
int report_errors()
|
||||
{
|
||||
if (detail::global_error_count > 0)
|
||||
@@ -277,4 +324,7 @@ int report_errors()
|
||||
#define CHECK_LE(X, Y) boost::math::test::check_le((X), (Y), __FILE__, __func__, __LINE__)
|
||||
|
||||
#define CHECK_CONDITIONED_ERROR(V, W, X, Y, Z) boost::math::test::check_conditioned_error((V), (W), (X), (Y), (Z), __FILE__, __func__, __LINE__)
|
||||
|
||||
#define CHECK_ABSOLUTE_ERROR(X, Y, Z) boost::math::test::check_absolute_error((X), (Y), (Z), __FILE__, __func__, __LINE__)
|
||||
|
||||
#endif
|
||||
|
||||
Reference in New Issue
Block a user