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)
7 #define BOOST_TEST_MODULE chebyshev_transform_test
9 #include <boost/cstdfloat.hpp>
10 #include <boost/type_index.hpp>
11 #include <boost/test/included/unit_test.hpp>
12 #include <boost/test/floating_point_comparison.hpp>
13 #include <boost/math/special_functions/chebyshev.hpp>
14 #include <boost/math/special_functions/chebyshev_transform.hpp>
15 #include <boost/math/special_functions/sinc.hpp>
16 #include <boost/multiprecision/cpp_bin_float.hpp>
17 #include <boost/multiprecision/cpp_dec_float.hpp>
19 #if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4)
26 using boost::multiprecision::cpp_bin_float_quad
;
27 using boost::multiprecision::cpp_bin_float_50
;
28 using boost::multiprecision::cpp_bin_float_100
;
29 using boost::math::chebyshev_t
;
30 using boost::math::chebyshev_t_prime
;
31 using boost::math::chebyshev_u
;
32 using boost::math::chebyshev_transform
;
36 void test_sin_chebyshev_transform()
38 using boost::math::chebyshev_transform
;
39 using boost::math::constants::half_pi
;
44 Real tol
= 10*std::numeric_limits
<Real
>::epsilon();
45 auto f
= [](Real x
) { return sin(x
); };
48 chebyshev_transform
<Real
> cheb(f
, a
, b
, tol
);
57 BOOST_CHECK_SMALL(cheb(x
), 100*tol
);
58 BOOST_CHECK_CLOSE_FRACTION(c
, cheb
.prime(x
), 100*tol
);
62 BOOST_CHECK_CLOSE_FRACTION(s
, cheb(x
), 100*tol
);
65 BOOST_CHECK_SMALL(cheb
.prime(x
), 100*tol
);
69 BOOST_CHECK_CLOSE_FRACTION(c
, cheb
.prime(x
), 100*tol
);
72 x
+= static_cast<Real
>(1)/static_cast<Real
>(1 << 7);
75 Real Q
= cheb
.integrate();
77 BOOST_CHECK_CLOSE_FRACTION(1 - cos(static_cast<Real
>(1)), Q
, 100*tol
);
82 void test_sinc_chebyshev_transform()
87 using boost::math::sinc_pi
;
88 using boost::math::chebyshev_transform
;
89 using boost::math::constants::half_pi
;
91 Real tol
= 500*std::numeric_limits
<Real
>::epsilon();
92 auto f
= [](Real x
) { return boost::math::sinc_pi(x
); };
95 chebyshev_transform
<Real
> cheb(f
, a
, b
, tol
/50);
101 Real ds
= (cos(x
)-sinc_pi(x
))/x
;
102 if (x
== 0) { ds
= 0; }
105 BOOST_CHECK_SMALL(cheb(x
), tol
);
109 BOOST_CHECK_CLOSE_FRACTION(s
, cheb(x
), tol
);
114 BOOST_CHECK_SMALL(cheb
.prime(x
), 5 * tol
);
118 BOOST_CHECK_CLOSE_FRACTION(ds
, cheb
.prime(x
), 300*tol
);
120 x
+= static_cast<Real
>(1)/static_cast<Real
>(1 << 7);
123 Real Q
= cheb
.integrate();
124 //NIntegrate[Sinc[x], {x, 0, 1}, WorkingPrecision -> 200, AccuracyGoal -> 150, PrecisionGoal -> 150, MaxRecursion -> 150]
125 Real Q_exp
= boost::lexical_cast
<Real
>("0.94608307036718301494135331382317965781233795473811179047145477356668");
126 BOOST_CHECK_CLOSE_FRACTION(Q_exp
, Q
, tol
);
131 //Examples taken from "Approximation Theory and Approximation Practice", by Trefethen
133 void test_atap_examples()
136 using boost::math::constants::half
;
137 using boost::math::sinc_pi
;
138 using boost::math::chebyshev_transform
;
139 using boost::math::constants::half_pi
;
141 Real tol
= 10*std::numeric_limits
<Real
>::epsilon();
142 auto f1
= [](Real x
) { return ((0 < x
) - (x
< 0)) - x
/2; };
143 auto f2
= [](Real x
) { Real t
= sin(6*x
); Real s
= sin(x
+ exp(2*x
));
144 Real u
= (0 < s
) - (s
< 0);
147 auto f3
= [](Real x
) { return sin(6*x
) + sin(60*exp(x
)); };
149 auto f4
= [](Real x
) { return 1/(1+1000*(x
+half
<Real
>())*(x
+half
<Real
>())) + 1/sqrt(1+1000*(x
-.5)*(x
-0.5));};
152 chebyshev_transform
<Real
> cheb1(f1
, a
, b
);
153 chebyshev_transform
<Real
> cheb2(f2
, a
, b
, tol
);
154 //chebyshev_transform<Real> cheb3(f3, a, b, tol);
160 if (sizeof(Real
) == sizeof(float))
162 BOOST_CHECK_CLOSE_FRACTION(f1(x
), cheb1(x
), 4e-3);
166 BOOST_CHECK_CLOSE_FRACTION(f1(x
), cheb1(x
), 1.3e-5);
168 BOOST_CHECK_CLOSE_FRACTION(f2(x
), cheb2(x
), 4e-3);
169 //BOOST_CHECK_CLOSE_FRACTION(f3(x), cheb3(x), 100*tol);
170 x
+= static_cast<Real
>(1)/static_cast<Real
>(1 << 7);
174 //Validate that the Chebyshev polynomials are well approximated by the Chebyshev transform.
176 void test_chebyshev_chebyshev_transform()
178 Real tol
= 500*std::numeric_limits
<Real
>::epsilon();
180 auto t0
= [](Real
) { return 1; };
181 chebyshev_transform
<Real
> cheb0(t0
, -1, 1);
182 BOOST_CHECK_CLOSE_FRACTION(cheb0
.coefficients()[0], 2, tol
);
187 BOOST_CHECK_CLOSE_FRACTION(cheb0(x
), 1, tol
);
188 BOOST_CHECK_SMALL(cheb0
.prime(x
), tol
);
189 x
+= static_cast<Real
>(1)/static_cast<Real
>(1 << 7);
193 auto t1
= [](Real x
) { return x
; };
194 chebyshev_transform
<Real
> cheb1(t1
, -1, 1);
195 BOOST_CHECK_CLOSE_FRACTION(cheb1
.coefficients()[1], 1, tol
);
202 BOOST_CHECK_SMALL(cheb1(x
), tol
);
206 BOOST_CHECK_CLOSE_FRACTION(cheb1(x
), x
, tol
);
208 BOOST_CHECK_CLOSE_FRACTION(cheb1
.prime(x
), 1, tol
);
209 x
+= static_cast<Real
>(1)/static_cast<Real
>(1 << 7);
213 auto t2
= [](Real x
) { return 2*x
*x
-1; };
214 chebyshev_transform
<Real
> cheb2(t2
, -1, 1);
215 BOOST_CHECK_CLOSE_FRACTION(cheb2
.coefficients()[2], 1, tol
);
220 BOOST_CHECK_CLOSE_FRACTION(cheb2(x
), t2(x
), tol
);
223 BOOST_CHECK_CLOSE_FRACTION(cheb2
.prime(x
), 4*x
, tol
);
227 BOOST_CHECK_SMALL(cheb2
.prime(x
), tol
);
229 x
+= static_cast<Real
>(1)/static_cast<Real
>(1 << 7);
233 BOOST_AUTO_TEST_CASE(chebyshev_transform_test
)
236 test_chebyshev_chebyshev_transform
<float>();
237 test_sin_chebyshev_transform
<float>();
238 test_atap_examples
<float>();
239 test_sinc_chebyshev_transform
<float>();
242 test_chebyshev_chebyshev_transform
<double>();
243 test_sin_chebyshev_transform
<double>();
244 test_atap_examples
<double>();
245 test_sinc_chebyshev_transform
<double>();
248 test_chebyshev_chebyshev_transform
<long double>();
249 test_sin_chebyshev_transform
<long double>();
250 test_atap_examples
<long double>();
251 test_sinc_chebyshev_transform
<long double>();
254 #ifdef BOOST_HAS_FLOAT128
255 test_chebyshev_chebyshev_transform
<__float128
>();
256 test_sin_chebyshev_transform
<__float128
>();
257 test_atap_examples
<__float128
>();
258 test_sinc_chebyshev_transform
<__float128
>();