libs/math/test/chebyshev_test.cpp

   1 /*
   2  * Copyright Nick Thompson, 2017
   3  * Use, modification and distribution are subject to the
   4  * Boost Software License, Version 1.0. (See accompanying file
   5  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6  */
   7
   8 #define BOOST_TEST_MODULE chebyshev_test
   9
  10 #include <boost/type_index.hpp>
  11 #include <boost/test/included/unit_test.hpp>
  12 #include <boost/test/tools/floating_point_comparison.hpp>
  13 #include <boost/math/special_functions/chebyshev.hpp>
  14 #include <boost/math/special_functions/sinc.hpp>
  15 #include <boost/multiprecision/cpp_bin_float.hpp>
  16 #include <boost/multiprecision/cpp_dec_float.hpp>
  17 #include <boost/array.hpp>
  18
  19 using boost::multiprecision::cpp_bin_float_quad;
  20 using boost::multiprecision::cpp_bin_float_50;
  21 using boost::multiprecision::cpp_bin_float_100;
  22 using boost::math::chebyshev_t;
  23 using boost::math::chebyshev_t_prime;
  24 using boost::math::chebyshev_u;
  25
  26 template<class Real>
  27 void test_polynomials()
  28 {
  29     std::cout << "Testing explicit polynomial representations of the Chebyshev polynomials on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
  30
  31     Real x = -2;
  32     Real tol = 400*std::numeric_limits<Real>::epsilon();
  33     if (tol > std::numeric_limits<float>::epsilon())
  34        tol *= 10;   // float results have much larger error rates.
  35     while (x < 2)
  36     {
  37         BOOST_CHECK_CLOSE_FRACTION(chebyshev_t(0, x), Real(1), tol);
  38         BOOST_CHECK_CLOSE_FRACTION(chebyshev_t(1, x), x, tol);
  39         BOOST_CHECK_CLOSE_FRACTION(chebyshev_t(2, x), 2*x*x - 1, tol);
  40         BOOST_CHECK_CLOSE_FRACTION(chebyshev_t(3, x), x*(4*x*x-3), tol);
  41         BOOST_CHECK_CLOSE_FRACTION(chebyshev_t(4, x), 8*x*x*(x*x - 1) + 1, tol);
  42         BOOST_CHECK_CLOSE_FRACTION(chebyshev_t(5, x), x*(16*x*x*x*x - 20*x*x + 5), tol);
  43         x += 1/static_cast<Real>(1<<7);
  44     }
  45
  46     x = -2;
  47     tol = 10*tol;
  48     while (x < 2)
  49     {
  50         BOOST_CHECK_CLOSE_FRACTION(chebyshev_u(0, x), Real(1), tol);
  51         BOOST_CHECK_CLOSE_FRACTION(chebyshev_u(1, x), 2*x, tol);
  52         BOOST_CHECK_CLOSE_FRACTION(chebyshev_u(2, x), 4*x*x - 1, tol);
  53         BOOST_CHECK_CLOSE_FRACTION(chebyshev_u(3, x), 4*x*(2*x*x - 1), tol);
  54         x += 1/static_cast<Real>(1<<7);
  55     }
  56 }
  57
  58
  59 template<class Real>
  60 void test_derivatives()
  61 {
  62     std::cout << "Testing explicit polynomial representations of the Chebyshev polynomial derivatives on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
  63
  64     Real x = -2;
  65     Real tol = 1000*std::numeric_limits<Real>::epsilon();
  66     while (x < 2)
  67     {
  68         BOOST_CHECK_CLOSE_FRACTION(chebyshev_t_prime(0, x), Real(0), tol);
  69         BOOST_CHECK_CLOSE_FRACTION(chebyshev_t_prime(1, x), Real(1), tol);
  70         BOOST_CHECK_CLOSE_FRACTION(chebyshev_t_prime(2, x), 4*x, tol);
  71         BOOST_CHECK_CLOSE_FRACTION(chebyshev_t_prime(3, x), 3*(4*x*x - 1), tol);
  72         BOOST_CHECK_CLOSE_FRACTION(chebyshev_t_prime(4, x), 16*x*(2*x*x - 1), tol);
  73         // This one makes the tolerance have to grow too large; the Chebyshev recurrence is more stable than naive polynomial evaluation anyway.
  74         //BOOST_CHECK_CLOSE_FRACTION(chebyshev_t_prime(5, x), 5*(4*x*x*(4*x*x - 3) + 1), tol);
  75         x += 1/static_cast<Real>(1<<7);
  76     }
  77 }
  78
  79 template<class Real>
  80 void test_clenshaw_recurrence()
  81 {
  82     using boost::math::chebyshev_clenshaw_recurrence;
  83     boost::array<Real, 5> c0 = { {2, 0, 0, 0, 0} };
  84     // Check the size = 1 case:
  85     boost::array<Real, 1> c01 = { {2} };
  86     // Check the size = 2 case:
  87     boost::array<Real, 2> c02 = { {2, 0} };
  88     boost::array<Real, 4> c1 = { {0, 1, 0, 0} };
  89     boost::array<Real, 4> c2 = { {0, 0, 1, 0} };
  90     boost::array<Real, 5> c3 = { {0, 0, 0, 1, 0} };
  91     boost::array<Real, 5> c4 = { {0, 0, 0, 0, 1} };
  92     boost::array<Real, 6> c5 = { {0, 0, 0, 0, 0, 1} };
  93     boost::array<Real, 7> c6 = { {0, 0, 0, 0, 0, 0, 1} };
  94
  95     Real x = -1;
  96     Real tol = 10*std::numeric_limits<Real>::epsilon();
  97     if (tol > std::numeric_limits<float>::epsilon())
  98        tol *= 100;   // float results have much larger error rates.
  99     while (x <= 1)
 100     {
 101         Real y = chebyshev_clenshaw_recurrence(c0.data(), c0.size(), x);
 102         BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(0, x), tol);
 103
 104         y = chebyshev_clenshaw_recurrence(c01.data(), c01.size(), x);
 105         BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(0, x), tol);
 106
 107         y = chebyshev_clenshaw_recurrence(c02.data(), c02.size(), x);
 108         BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(0, x), tol);
 109
 110         y = chebyshev_clenshaw_recurrence(c1.data(), c1.size(), x);
 111         BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(1, x), tol);
 112
 113         y = chebyshev_clenshaw_recurrence(c2.data(), c2.size(), x);
 114         BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(2, x), tol);
 115
 116         y = chebyshev_clenshaw_recurrence(c3.data(), c3.size(), x);
 117         BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(3, x), tol);
 118
 119         y = chebyshev_clenshaw_recurrence(c4.data(), c4.size(), x);
 120         BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(4, x), tol);
 121
 122         y = chebyshev_clenshaw_recurrence(c5.data(), c5.size(), x);
 123         BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(5, x), tol);
 124
 125         y = chebyshev_clenshaw_recurrence(c6.data(), c6.size(), x);
 126         BOOST_CHECK_CLOSE_FRACTION(y, chebyshev_t(6, x), tol);
 127
 128         x += static_cast<Real>(1)/static_cast<Real>(1 << 7);
 129     }
 130 }
 131
 132 BOOST_AUTO_TEST_CASE(chebyshev_test)
 133 {
 134     test_clenshaw_recurrence<float>();
 135     test_clenshaw_recurrence<double>();
 136     test_clenshaw_recurrence<long double>();
 137
 138     test_polynomials<float>();
 139     test_polynomials<double>();
 140     test_polynomials<long double>();
 141     test_polynomials<cpp_bin_float_quad>();
 142
 143     test_derivatives<float>();
 144     test_derivatives<double>();
 145     test_derivatives<long double>();
 146     test_derivatives<cpp_bin_float_quad>();
 147 }