2 * Copyright Nick Thompson, 2019
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)
8 #include "math_unit_test.hpp"
11 #include <boost/math/constants/constants.hpp>
12 #include <boost/math/interpolators/cardinal_trigonometric.hpp>
13 #ifdef BOOST_HAS_FLOAT128
14 #include <boost/multiprecision/float128.hpp>
19 using boost::math::constants::two_pi
;
20 using boost::math::interpolators::cardinal_trigonometric
;
27 for(size_t n
= 1; n
< 20; ++n
)
30 std::vector
<Real
> v(n
, c
);
31 auto ct
= cardinal_trigonometric
<decltype(v
)>(v
, t0
, h
);
32 CHECK_ULP_CLOSE(c
, ct(0.3), 3);
33 CHECK_ULP_CLOSE(c
*h
*n
, ct
.integrate(), 3);
34 CHECK_ULP_CLOSE(c
*c
*h
*n
, ct
.squared_l2(), 3);
35 CHECK_MOLLIFIED_CLOSE(Real(0), ct
.prime(0.8), 25*std::numeric_limits
<Real
>::epsilon());
36 CHECK_MOLLIFIED_CLOSE(Real(0), ct
.double_prime(0.8), 25*std::numeric_limits
<Real
>::epsilon());
41 void test_interpolation_condition()
43 std::mt19937
gen(1234);
44 std::uniform_real_distribution
<Real
> dis(1, 10);
46 for(size_t n
= 1; n
< 20; ++n
) {
49 std::vector
<Real
> v(n
);
50 for (size_t i
= 0; i
< n
; ++i
) {
53 auto ct
= cardinal_trigonometric
<decltype(v
)>(v
, t0
, h
);
54 for (size_t i
= 0; i
< n
; ++i
) {
57 Real computed
= ct(arg
);
58 if(!CHECK_ULP_CLOSE(expected
, computed
, 5*n
))
60 std::cerr
<< " Samples: " << n
<< "\n";
68 #ifdef BOOST_HAS_FLOAT128
69 void test_constant_q()
73 for(size_t n
= 1; n
< 20; ++n
)
76 std::vector
<__float128
> v(n
, c
);
77 auto ct
= cardinal_trigonometric
<decltype(v
)>(v
, t0
, h
);
78 CHECK_ULP_CLOSE(boost::multiprecision::float128(c
), boost::multiprecision::float128(ct(0.3)), 3);
79 CHECK_ULP_CLOSE(boost::multiprecision::float128(c
*h
*n
), boost::multiprecision::float128(ct
.integrate()), 3);
86 void test_sampled_sine()
90 for (unsigned n
= 15; n
< 50; ++n
)
95 std::vector
<Real
> v(n
);
96 auto s
= [&](Real t
) { return sin(two_pi
<Real
>()*(t
-t0
)/T
);};
97 auto s_prime
= [&](Real t
) { return two_pi
<Real
>()*cos(two_pi
<Real
>()*(t
-t0
)/T
)/T
;};
98 auto s_double_prime
= [&](Real t
) { return -two_pi
<Real
>()*two_pi
<Real
>()*sin(two_pi
<Real
>()*(t
-t0
)/T
)/(T
*T
);};
99 for(size_t j
= 0; j
< v
.size(); ++j
)
104 auto ct
= cardinal_trigonometric
<decltype(v
)>(v
, t0
, h
);
105 CHECK_ULP_CLOSE(T
, ct
.period(), 3);
106 std::mt19937
gen(1234);
107 std::uniform_real_distribution
<Real
> dist(0, 500);
111 Real arg
= dist(gen
);
112 Real expected
= s(arg
);
113 Real computed
= ct(arg
);
114 CHECK_MOLLIFIED_CLOSE(expected
, computed
, std::numeric_limits
<Real
>::epsilon()*4000);
116 expected
= s_prime(arg
);
117 computed
= ct
.prime(arg
);
118 CHECK_MOLLIFIED_CLOSE(expected
, computed
, 18000*std::numeric_limits
<Real
>::epsilon());
120 expected
= s_double_prime(arg
);
121 computed
= ct
.double_prime(arg
);
122 CHECK_MOLLIFIED_CLOSE(expected
, computed
, 100000*std::numeric_limits
<Real
>::epsilon());
125 CHECK_MOLLIFIED_CLOSE(Real(0), ct
.integrate(), std::numeric_limits
<Real
>::epsilon());
136 auto bump
= [](Real x
)->Real
{ if (abs(x
) >= 1) { return Real(0); } return exp(-Real(1)/(Real(1)-x
*x
)); };
137 auto bump_prime
= [](Real x
)->Real
{
138 if (abs(x
) >= 1) { return Real(0); }
140 return -2*x
*exp(-Real(1)/(Real(1)-x
*x
))/pow(1-x
*x
,2);
143 auto bump_double_prime
= [](Real x
)->Real
{
144 if (abs(x
) >= 1) { return Real(0); }
146 return (6*pow(x
,4)-2)*exp(-Real(1)/(Real(1)-x
*x
))/pow(1-x
*x
,4);
152 Real h
= Real(2)/Real(n
);
154 std::vector
<Real
> v(n
);
155 for(size_t i
= 0; i
< n
; ++i
)
161 auto ct
= cardinal_trigonometric
<decltype(v
)>(v
, t0
, h
);
162 std::mt19937
gen(323723);
163 std::uniform_real_distribution
<long double> dis(-0.9, 0.9);
168 Real t
= static_cast<Real
>(dis(gen
));
169 Real expected
= bump(t
);
170 Real computed
= ct(t
);
171 if(!CHECK_MOLLIFIED_CLOSE(expected
, computed
, 2*std::numeric_limits
<Real
>::epsilon())) {
172 std::cerr
<< " Problem occurred at abscissa " << t
<< "\n";
175 expected
= bump_prime(t
);
176 computed
= ct
.prime(t
);
177 if(!CHECK_MOLLIFIED_CLOSE(expected
, computed
, 4000*std::numeric_limits
<Real
>::epsilon())) {
178 std::cerr
<< " Problem occurred at abscissa " << t
<< "\n";
181 expected
= bump_double_prime(t
);
182 computed
= ct
.double_prime(t
);
183 if(!CHECK_MOLLIFIED_CLOSE(expected
, computed
, 4000*4000*std::numeric_limits
<Real
>::epsilon())) {
184 std::cerr
<< " Problem occurred at abscissa " << t
<< "\n";
191 // NIntegrate[Exp[-1/(1-x*x)],{x,-1,1}]
192 CHECK_ULP_CLOSE(Real(0.443993816168079437823L), ct
.integrate(), 3);
194 // NIntegrate[Exp[-2/(1-x*x)],{x,-1,1}]
195 CHECK_ULP_CLOSE(Real(0.1330861208449942715569473279553285713625791551628130055345002588895389L), ct
.squared_l2(), 1);
205 test_constant
<float>();
206 test_sampled_sine
<float>();
208 test_interpolation_condition
<float>();
213 test_constant
<double>();
214 test_sampled_sine
<double>();
216 test_interpolation_condition
<double>();
220 test_constant
<long double>();
221 test_sampled_sine
<long double>();
222 test_bump
<long double>();
223 test_interpolation_condition
<long double>();
227 #ifdef BOOST_HAS_FLOAT128
232 return boost::math::test::report_errors();