1 /*
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)
6 */
7
8 #include "math_unit_test.hpp"
9 #include <numeric>
10 #include <utility>
11 #include <boost/math/interpolators/cardinal_quadratic_b_spline.hpp>
12 using boost::math::interpolators::cardinal_quadratic_b_spline;
13
14 template<class Real>
test_constant()15 void test_constant()
16 {
17 Real c = 7.2;
18 Real t0 = 0;
19 Real h = Real(1)/Real(16);
20 size_t n = 512;
21 std::vector<Real> v(n, c);
22 auto qbs = cardinal_quadratic_b_spline<Real>(v.data(), v.size(), t0, h);
23
24 size_t i = 0;
25 while (i < n) {
26 Real t = t0 + i*h;
27 CHECK_ULP_CLOSE(c, qbs(t), 2);
28 CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 100*std::numeric_limits<Real>::epsilon());
29 ++i;
30 }
31
32 i = 0;
33 while (i < n) {
34 Real t = t0 + i*h + h/2;
35 CHECK_ULP_CLOSE(c, qbs(t), 2);
36 CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 300*std::numeric_limits<Real>::epsilon());
37 t = t0 + i*h + h/4;
38 CHECK_ULP_CLOSE(c, qbs(t), 2);
39 CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 150*std::numeric_limits<Real>::epsilon());
40 ++i;
41 }
42 }
43
44 template<class Real>
test_linear()45 void test_linear()
46 {
47 Real m = 8.3;
48 Real b = 7.2;
49 Real t0 = 0;
50 Real h = Real(1)/Real(16);
51 size_t n = 512;
52 std::vector<Real> y(n);
53 for (size_t i = 0; i < n; ++i) {
54 Real t = i*h;
55 y[i] = m*t + b;
56 }
57 auto qbs = cardinal_quadratic_b_spline<Real>(y.data(), y.size(), t0, h);
58
59 size_t i = 0;
60 while (i < n) {
61 Real t = t0 + i*h;
62 CHECK_ULP_CLOSE(m*t+b, qbs(t), 2);
63 CHECK_ULP_CLOSE(m, qbs.prime(t), 820);
64 ++i;
65 }
66
67 i = 0;
68 while (i < n) {
69 Real t = t0 + i*h + h/2;
70 CHECK_ULP_CLOSE(m*t+b, qbs(t), 2);
71 CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits<Real>::epsilon());
72 t = t0 + i*h + h/4;
73 CHECK_ULP_CLOSE(m*t+b, qbs(t), 3);
74 CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits<Real>::epsilon());
75 ++i;
76 }
77 }
78
79 template<class Real>
test_quadratic()80 void test_quadratic()
81 {
82 Real a = 8.2;
83 Real b = 7.2;
84 Real c = -9.2;
85 Real t0 = 0;
86 Real h = Real(1)/Real(16);
87 size_t n = 513;
88 std::vector<Real> y(n);
89 for (size_t i = 0; i < n; ++i) {
90 Real t = i*h;
91 y[i] = a*t*t + b*t + c;
92 }
93 Real t_max = t0 + (n-1)*h;
94 auto qbs = cardinal_quadratic_b_spline<Real>(y, t0, h, b, 2*a*t_max + b);
95
96 size_t i = 0;
97 while (i < n) {
98 Real t = t0 + i*h;
99 CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 2);
100 ++i;
101 }
102
103 i = 0;
104 while (i < n) {
105 Real t = t0 + i*h + h/2;
106 CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 47);
107
108 t = t0 + i*h + h/4;
109 if (!CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 104)) {
110 std::cerr << " Problem abscissa t = " << t << "\n";
111 }
112 ++i;
113 }
114 }
115
main()116 int main()
117 {
118 test_constant<float>();
119 test_constant<double>();
120 test_constant<long double>();
121
122 test_linear<float>();
123 test_linear<double>();
124 test_linear<long double>();
125
126 test_quadratic<double>();
127 test_quadratic<long double>();
128
129 return boost::math::test::report_errors();
130 }
131