]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
1 | // (C) Copyright Nick Thompson 2019. |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | ||
6 | #ifndef BOOST_MATH_SPECIAL_CARDINAL_B_SPLINE_HPP | |
7 | #define BOOST_MATH_SPECIAL_CARDINAL_B_SPLINE_HPP | |
8 | ||
9 | #include <array> | |
10 | #include <cmath> | |
11 | #include <limits> | |
12 | #include <type_traits> | |
13 | ||
14 | namespace boost { namespace math { | |
15 | ||
16 | namespace detail { | |
17 | ||
18 | template<class Real> | |
19 | inline Real B1(Real x) | |
20 | { | |
21 | if (x < 0) | |
22 | { | |
23 | return B1(-x); | |
24 | } | |
25 | if (x < Real(1)) | |
26 | { | |
27 | return 1 - x; | |
28 | } | |
29 | return Real(0); | |
30 | } | |
31 | } | |
32 | ||
33 | template<unsigned n, typename Real> | |
34 | Real cardinal_b_spline(Real x) { | |
35 | static_assert(!std::is_integral<Real>::value, "Does not work with integral types."); | |
36 | ||
37 | if (x < 0) { | |
38 | // All B-splines are even functions: | |
39 | return cardinal_b_spline<n, Real>(-x); | |
40 | } | |
41 | ||
42 | if (n==0) | |
43 | { | |
44 | if (x < Real(1)/Real(2)) { | |
45 | return Real(1); | |
46 | } | |
47 | else if (x == Real(1)/Real(2)) { | |
48 | return Real(1)/Real(2); | |
49 | } | |
50 | else { | |
51 | return Real(0); | |
52 | } | |
53 | } | |
54 | ||
55 | if (n==1) | |
56 | { | |
57 | return detail::B1(x); | |
58 | } | |
59 | ||
60 | Real supp_max = (n+1)/Real(2); | |
61 | if (x >= supp_max) | |
62 | { | |
63 | return Real(0); | |
64 | } | |
65 | ||
66 | // Fill v with values of B1: | |
67 | // At most two of these terms are nonzero, and at least 1. | |
68 | // There is only one non-zero term when n is odd and x = 0. | |
69 | std::array<Real, n> v; | |
70 | Real z = x + 1 - supp_max; | |
71 | for (unsigned i = 0; i < n; ++i) | |
72 | { | |
73 | v[i] = detail::B1(z); | |
74 | z += 1; | |
75 | } | |
76 | ||
77 | Real smx = supp_max - x; | |
78 | for (unsigned j = 2; j <= n; ++j) | |
79 | { | |
80 | Real a = (j + 1 - smx); | |
81 | Real b = smx; | |
82 | for(unsigned k = 0; k <= n - j; ++k) | |
83 | { | |
84 | v[k] = (a*v[k+1] + b*v[k])/Real(j); | |
85 | a += 1; | |
86 | b -= 1; | |
87 | } | |
88 | } | |
89 | ||
90 | return v[0]; | |
91 | } | |
92 | ||
93 | ||
94 | template<unsigned n, typename Real> | |
95 | Real cardinal_b_spline_prime(Real x) | |
96 | { | |
97 | static_assert(!std::is_integral<Real>::value, "Cardinal B-splines do not work with integer types."); | |
98 | ||
99 | if (x < 0) | |
100 | { | |
101 | // All B-splines are even functions, so derivatives are odd: | |
102 | return -cardinal_b_spline_prime<n, Real>(-x); | |
103 | } | |
104 | ||
105 | ||
106 | if (n==0) | |
107 | { | |
108 | // Kinda crazy but you get what you ask for! | |
109 | if (x == Real(1)/Real(2)) | |
110 | { | |
111 | return std::numeric_limits<Real>::infinity(); | |
112 | } | |
113 | else | |
114 | { | |
115 | return Real(0); | |
116 | } | |
117 | } | |
118 | ||
119 | if (n==1) | |
120 | { | |
121 | if (x==0) | |
122 | { | |
123 | return Real(0); | |
124 | } | |
125 | if (x==1) | |
126 | { | |
127 | return -Real(1)/Real(2); | |
128 | } | |
129 | return Real(-1); | |
130 | } | |
131 | ||
132 | ||
133 | Real supp_max = (n+1)/Real(2); | |
134 | if (x >= supp_max) | |
135 | { | |
136 | return Real(0); | |
137 | } | |
138 | ||
139 | // Now we want to evaluate B_{n}(x), but stop at the second to last step and collect B_{n-1}(x+1/2) and B_{n-1}(x-1/2): | |
140 | std::array<Real, n> v; | |
141 | Real z = x + 1 - supp_max; | |
142 | for (unsigned i = 0; i < n; ++i) | |
143 | { | |
144 | v[i] = detail::B1(z); | |
145 | z += 1; | |
146 | } | |
147 | ||
148 | Real smx = supp_max - x; | |
149 | for (unsigned j = 2; j <= n - 1; ++j) | |
150 | { | |
151 | Real a = (j + 1 - smx); | |
152 | Real b = smx; | |
153 | for(unsigned k = 0; k <= n - j; ++k) | |
154 | { | |
155 | v[k] = (a*v[k+1] + b*v[k])/Real(j); | |
156 | a += 1; | |
157 | b -= 1; | |
158 | } | |
159 | } | |
160 | ||
161 | return v[1] - v[0]; | |
162 | } | |
163 | ||
164 | ||
165 | template<unsigned n, typename Real> | |
166 | Real cardinal_b_spline_double_prime(Real x) | |
167 | { | |
168 | static_assert(!std::is_integral<Real>::value, "Cardinal B-splines do not work with integer types."); | |
169 | static_assert(n >= 3, "n>=3 for second derivatives of cardinal B-splines is required."); | |
170 | ||
171 | if (x < 0) | |
172 | { | |
173 | // All B-splines are even functions, so second derivatives are even: | |
174 | return cardinal_b_spline_double_prime<n, Real>(-x); | |
175 | } | |
176 | ||
177 | ||
178 | Real supp_max = (n+1)/Real(2); | |
179 | if (x >= supp_max) | |
180 | { | |
181 | return Real(0); | |
182 | } | |
183 | ||
184 | // Now we want to evaluate B_{n}(x), but stop at the second to last step and collect B_{n-1}(x+1/2) and B_{n-1}(x-1/2): | |
185 | std::array<Real, n> v; | |
186 | Real z = x + 1 - supp_max; | |
187 | for (unsigned i = 0; i < n; ++i) | |
188 | { | |
189 | v[i] = detail::B1(z); | |
190 | z += 1; | |
191 | } | |
192 | ||
193 | Real smx = supp_max - x; | |
194 | for (unsigned j = 2; j <= n - 2; ++j) | |
195 | { | |
196 | Real a = (j + 1 - smx); | |
197 | Real b = smx; | |
198 | for(unsigned k = 0; k <= n - j; ++k) | |
199 | { | |
200 | v[k] = (a*v[k+1] + b*v[k])/Real(j); | |
201 | a += 1; | |
202 | b -= 1; | |
203 | } | |
204 | } | |
205 | ||
206 | return v[2] - 2*v[1] + v[0]; | |
207 | } | |
208 | ||
209 | ||
210 | template<unsigned n, class Real> | |
211 | Real forward_cardinal_b_spline(Real x) | |
212 | { | |
213 | static_assert(!std::is_integral<Real>::value, "Cardinal B-splines do not work with integral types."); | |
214 | return cardinal_b_spline<n>(x - (n+1)/Real(2)); | |
215 | } | |
216 | ||
217 | }} | |
218 | #endif |