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"
13 #include <boost/math/special_functions/cardinal_b_spline.hpp>
14 #include <boost/math/interpolators/detail/cubic_b_spline_detail.hpp>
15 #ifdef BOOST_HAS_FLOAT128
16 #include <boost/multiprecision/float128.hpp>
17 using boost::multiprecision::float128
;
21 using boost::math::cardinal_b_spline
;
22 using boost::math::cardinal_b_spline_prime
;
23 using boost::math::forward_cardinal_b_spline
;
24 using boost::math::cardinal_b_spline_double_prime
;
30 Real t
= cardinal_b_spline
<0>(Real(1.1));
32 CHECK_ULP_CLOSE(expected
, t
, 0);
33 CHECK_ULP_CLOSE(expected
, cardinal_b_spline_prime
<0>(Real(1.1)), 0);
35 t
= cardinal_b_spline
<0>(Real(-1.1));
37 CHECK_ULP_CLOSE(expected
, t
, 0);
39 Real h
= Real(1)/Real(256);
40 for (t
= -Real(1)/Real(2)+h
; t
< Real(1)/Real(2); t
+= h
)
43 CHECK_ULP_CLOSE(expected
, cardinal_b_spline
<0>(t
), 0);
45 CHECK_ULP_CLOSE(expected
, cardinal_b_spline_prime
<0>(Real(1.1)), 0);
48 for (t
= h
; t
< 1; t
+= h
)
51 CHECK_ULP_CLOSE(expected
, forward_cardinal_b_spline
<0>(t
), 0);
58 Real t
= cardinal_b_spline
<1>(Real(2.1));
60 CHECK_ULP_CLOSE(expected
, t
, 0);
62 t
= cardinal_b_spline
<1>(Real(-2.1));
64 CHECK_ULP_CLOSE(expected
, t
, 0);
66 Real h
= Real(1)/Real(256);
67 for (t
= -1; t
<= 1; t
+= h
)
70 if(!CHECK_ULP_CLOSE(expected
, cardinal_b_spline
<1>(t
), 0) )
72 std::cerr
<< " Problem at t = " << t
<< "\n";
75 if (!CHECK_ULP_CLOSE(Real(1)/Real(2), cardinal_b_spline_prime
<1>(t
), 0)) {
76 std::cout
<< " Problem at t = " << t
<< "\n";
79 else if (t
== Real(1)) {
80 CHECK_ULP_CLOSE(-Real(1)/Real(2), cardinal_b_spline_prime
<1>(t
), 0);
83 CHECK_ULP_CLOSE(Real(1), cardinal_b_spline_prime
<1>(t
), 0);
86 CHECK_ULP_CLOSE(Real(0), cardinal_b_spline_prime
<1>(t
), 0);
89 CHECK_ULP_CLOSE(Real(-1), cardinal_b_spline_prime
<1>(t
), 0);
93 for (t
= 0; t
< 2; t
+= h
)
95 expected
= 1 - abs(t
-1);
96 CHECK_ULP_CLOSE(expected
, forward_cardinal_b_spline
<1>(t
), 0);
101 void test_quadratic()
104 auto b2
= [](Real x
) {
106 if (absx
>= 3/Real(2)) {
109 if (absx
>= 1/Real(2)) {
110 Real t
= absx
- 3/Real(2);
113 Real t1
= absx
- 1/Real(2);
114 Real t2
= absx
+ 1/Real(2);
115 return (2-t1
*t1
-t2
*t2
)/2;
118 auto b2_prime
= [&](Real x
)->Real
{
124 if (absx
>= 3/Real(2)) {
127 if (absx
>= 1/Real(2)) {
128 return (absx
- 3/Real(2))*signx
;
130 return -2*absx
*signx
;
134 Real h
= 1/Real(256);
135 for (Real t
= -5; t
<= 5; t
+= h
) {
136 Real expected
= b2(t
);
137 CHECK_ULP_CLOSE(expected
, cardinal_b_spline
<2>(t
), 0);
138 expected
= b2_prime(t
);
140 if (!CHECK_ULP_CLOSE(expected
, cardinal_b_spline_prime
<2>(t
), 0))
142 std::cerr
<< " Problem at t = " << t
<< "\n";
151 Real expected
= Real(2)/Real(3);
152 Real computed
= cardinal_b_spline
<3, Real
>(0);
153 CHECK_ULP_CLOSE(expected
, computed
, 0);
155 expected
= Real(1)/Real(6);
156 computed
= cardinal_b_spline
<3, Real
>(1);
157 CHECK_ULP_CLOSE(expected
, computed
, 0);
160 computed
= cardinal_b_spline
<3, Real
>(2);
161 CHECK_ULP_CLOSE(expected
, computed
, 0);
163 Real h
= 1/Real(256);
164 for (Real t
= -4; t
<= 4; t
+= h
) {
165 expected
= boost::math::detail::b3_spline_prime
<Real
>(t
);
166 computed
= cardinal_b_spline_prime
<3>(t
);
167 CHECK_ULP_CLOSE(expected
, computed
, 0);
168 expected
= boost::math::detail::b3_spline_double_prime
<Real
>(t
);
169 computed
= cardinal_b_spline_double_prime
<3>(t
);
170 if (!CHECK_ULP_CLOSE(expected
, computed
, 0)) {
171 std::cerr
<< " Problem at t = " << t
<< "\n";
179 Real expected
= Real(11)/Real(20);
180 Real computed
= cardinal_b_spline
<5, Real
>(0);
181 CHECK_ULP_CLOSE(expected
, computed
, 0);
183 expected
= Real(13)/Real(60);
184 computed
= cardinal_b_spline
<5, Real
>(1);
185 CHECK_ULP_CLOSE(expected
, computed
, 1);
187 expected
= Real(1)/Real(120);
188 computed
= cardinal_b_spline
<5, Real
>(2);
189 CHECK_ULP_CLOSE(expected
, computed
, 0);
192 computed
= cardinal_b_spline
<5, Real
>(3);
193 CHECK_ULP_CLOSE(expected
, computed
, 0);
197 template<unsigned n
, typename Real
>
198 void test_b_spline_derivatives()
200 Real h
= 1/Real(256);
201 Real supp
= (n
+Real(1))/Real(2);
202 for (Real t
= -supp
- 1; t
<= supp
+1; t
+= h
)
204 Real expected
= cardinal_b_spline
<n
-1>(t
+Real(1)/Real(2)) - cardinal_b_spline
<n
-1>(t
- Real(1)/Real(2));
205 Real computed
= cardinal_b_spline_prime
<n
>(t
);
206 CHECK_MOLLIFIED_CLOSE(expected
, computed
, std::numeric_limits
<Real
>::epsilon());
208 expected
= cardinal_b_spline
<n
-2>(t
+1) - 2*cardinal_b_spline
<n
-2>(t
) + cardinal_b_spline
<n
-2>(t
-1);
209 computed
= cardinal_b_spline_double_prime
<n
>(t
);
210 CHECK_MOLLIFIED_CLOSE(expected
, computed
, 2*std::numeric_limits
<Real
>::epsilon());
214 template<unsigned n
, typename Real
>
215 void test_partition_of_unity()
217 std::mt19937
gen(323723);
218 Real supp
= (n
+1.0)/2.0;
219 std::uniform_real_distribution
<Real
> dis(-supp
, -supp
+1);
221 for(size_t i
= 0; i
< 500; ++i
) {
225 one
+= cardinal_b_spline
<n
>(x
);
228 if(!CHECK_ULP_CLOSE(Real(1), one
, n
)) {
229 std::cerr
<< " Partition of unity failure at n = " << n
<< "\n";
239 test_box
<long double>();
243 test_hat
<long double>();
245 test_quadratic
<float>();
246 test_quadratic
<double>();
247 test_quadratic
<long double>();
250 test_cubic
<double>();
251 test_cubic
<long double>();
253 test_quintic
<float>();
254 test_quintic
<double>();
255 test_quintic
<long double>();
257 test_partition_of_unity
<1, double>();
258 test_partition_of_unity
<2, double>();
259 test_partition_of_unity
<3, double>();
260 test_partition_of_unity
<4, double>();
261 test_partition_of_unity
<5, double>();
262 test_partition_of_unity
<6, double>();
264 test_b_spline_derivatives
<3, double>();
265 test_b_spline_derivatives
<4, double>();
266 test_b_spline_derivatives
<5, double>();
267 test_b_spline_derivatives
<6, double>();
268 test_b_spline_derivatives
<7, double>();
269 test_b_spline_derivatives
<8, double>();
270 test_b_spline_derivatives
<9, double>();
272 test_b_spline_derivatives
<3, long double>();
273 test_b_spline_derivatives
<4, long double>();
274 test_b_spline_derivatives
<5, long double>();
275 test_b_spline_derivatives
<6, long double>();
276 test_b_spline_derivatives
<7, long double>();
277 test_b_spline_derivatives
<8, long double>();
278 test_b_spline_derivatives
<9, long double>();
281 #ifdef BOOST_HAS_FLOAT128
282 test_box
<float128
>();
283 test_hat
<float128
>();
284 test_quadratic
<float128
>();
285 test_cubic
<float128
>();
286 test_quintic
<float128
>();
289 return boost::math::test::report_errors();