]>
Commit | Line | Data |
---|---|---|
f67539c2 TL |
1 | /* |
2 | * Copyright Nick Thompson, 2020 | |
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 <vector> | |
12 | #include <array> | |
13 | #include <boost/random/uniform_real.hpp> | |
14 | #include <boost/random/mersenne_twister.hpp> | |
15 | #include <boost/math/interpolators/quintic_hermite.hpp> | |
16 | #include <boost/math/special_functions/next.hpp> | |
17 | #include <boost/circular_buffer.hpp> | |
18 | #ifdef BOOST_HAS_FLOAT128 | |
19 | #include <boost/multiprecision/float128.hpp> | |
20 | using boost::multiprecision::float128; | |
21 | #endif | |
22 | ||
23 | ||
24 | using boost::math::interpolators::quintic_hermite; | |
25 | using boost::math::interpolators::cardinal_quintic_hermite; | |
26 | using boost::math::interpolators::cardinal_quintic_hermite_aos; | |
27 | ||
28 | template<typename Real> | |
29 | void test_constant() | |
30 | { | |
31 | std::vector<Real> x{0,1,2,3, 9, 22, 81}; | |
32 | std::vector<Real> y(x.size()); | |
33 | std::vector<Real> dydx(x.size(), 0); | |
34 | std::vector<Real> d2ydx2(x.size(), 0); | |
35 | for (auto & t : y) | |
36 | { | |
37 | t = 7; | |
38 | } | |
39 | ||
40 | auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)); | |
41 | for (Real t = 0; t <= 81; t += 0.25) | |
42 | { | |
43 | CHECK_ULP_CLOSE(Real(7), qh(t), 24); | |
44 | CHECK_ULP_CLOSE(Real(0), qh.prime(t), 24); | |
45 | CHECK_ULP_CLOSE(Real(0), qh.double_prime(t), 24); | |
46 | } | |
47 | } | |
48 | ||
49 | ||
50 | template<typename Real> | |
51 | void test_linear() | |
52 | { | |
53 | std::vector<Real> x{0,1,2,3, 4,5,6,7,8,9}; | |
54 | std::vector<Real> y = x; | |
55 | std::vector<Real> dydx(x.size(), 1); | |
56 | std::vector<Real> d2ydx2(x.size(), 0); | |
57 | ||
58 | auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)); | |
59 | ||
60 | for (Real t = 0; t <= 9; t += 0.25) | |
61 | { | |
62 | CHECK_ULP_CLOSE(Real(t), qh(t), 2); | |
63 | CHECK_ULP_CLOSE(Real(1), qh.prime(t), 2); | |
64 | CHECK_ULP_CLOSE(Real(0), qh.double_prime(t), 2); | |
65 | } | |
66 | ||
67 | boost::random::mt19937 rng; | |
68 | boost::random::uniform_real_distribution<Real> dis(0.5,1); | |
69 | x.resize(512); | |
70 | x[0] = dis(rng); | |
71 | Real xmin = x[0]; | |
72 | for (size_t i = 1; i < x.size(); ++i) | |
73 | { | |
74 | x[i] = x[i-1] + dis(rng); | |
75 | } | |
76 | Real xmax = x.back(); | |
77 | ||
78 | y = x; | |
79 | dydx.resize(x.size(), 1); | |
80 | d2ydx2.resize(x.size(), 0); | |
81 | ||
82 | qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)); | |
83 | ||
84 | for (Real t = xmin; t <= xmax; t += 0.125) | |
85 | { | |
86 | CHECK_ULP_CLOSE(t, qh(t), 2); | |
87 | CHECK_ULP_CLOSE(Real(1), qh.prime(t), 100); | |
88 | CHECK_MOLLIFIED_CLOSE(Real(0), qh.double_prime(t), 200*std::numeric_limits<Real>::epsilon()); | |
89 | } | |
90 | } | |
91 | ||
92 | template<typename Real> | |
93 | void test_quadratic() | |
94 | { | |
95 | ||
96 | std::vector<Real> x{0,1,2,3, 4,5,6,7,8,9}; | |
97 | std::vector<Real> y(x.size()); | |
98 | for (size_t i = 0; i < y.size(); ++i) | |
99 | { | |
100 | y[i] = x[i]*x[i]/2; | |
101 | } | |
102 | ||
103 | std::vector<Real> dydx(x.size()); | |
104 | for (size_t i = 0; i < y.size(); ++i) { | |
105 | dydx[i] = x[i]; | |
106 | } | |
107 | ||
108 | std::vector<Real> d2ydx2(x.size(), 1); | |
109 | ||
110 | auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)); | |
111 | ||
112 | for (Real t = 0; t <= 9; t += 0.0078125) | |
113 | { | |
114 | CHECK_ULP_CLOSE(Real(t*t)/2, qh(t), 2); | |
115 | CHECK_ULP_CLOSE(t, qh.prime(t), 12); | |
116 | CHECK_ULP_CLOSE(Real(1), qh.double_prime(t), 32); | |
117 | } | |
118 | ||
119 | boost::random::mt19937 rng; | |
120 | boost::random::uniform_real_distribution<Real> dis(0.5,1); | |
121 | x.resize(8); | |
122 | x[0] = dis(rng); | |
123 | Real xmin = x[0]; | |
124 | for (size_t i = 1; i < x.size(); ++i) | |
125 | { | |
126 | x[i] = x[i-1] + dis(rng); | |
127 | } | |
128 | Real xmax = x.back(); | |
129 | ||
130 | y.resize(x.size()); | |
131 | for (size_t i = 0; i < y.size(); ++i) | |
132 | { | |
133 | y[i] = x[i]*x[i]/2; | |
134 | } | |
135 | ||
136 | dydx.resize(x.size()); | |
137 | for (size_t i = 0; i < y.size(); ++i) | |
138 | { | |
139 | dydx[i] = x[i]; | |
140 | } | |
141 | ||
142 | d2ydx2.resize(x.size(), 1); | |
143 | ||
144 | qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)); | |
145 | ||
146 | for (Real t = xmin; t <= xmax; t += 0.125) | |
147 | { | |
148 | CHECK_ULP_CLOSE(Real(t*t)/2, qh(t), 4); | |
149 | CHECK_ULP_CLOSE(t, qh.prime(t), 53); | |
150 | CHECK_ULP_CLOSE(Real(1), qh.double_prime(t), 700); | |
151 | } | |
152 | } | |
153 | ||
154 | template<typename Real> | |
155 | void test_cubic() | |
156 | { | |
157 | ||
158 | std::vector<Real> x{0,1,2,3, 4,5,6,7,8,9}; | |
159 | std::vector<Real> y(x.size()); | |
160 | for (size_t i = 0; i < y.size(); ++i) | |
161 | { | |
162 | y[i] = x[i]*x[i]*x[i]; | |
163 | } | |
164 | ||
165 | std::vector<Real> dydx(x.size()); | |
166 | for (size_t i = 0; i < y.size(); ++i) { | |
167 | dydx[i] = 3*x[i]*x[i]; | |
168 | } | |
169 | ||
170 | std::vector<Real> d2ydx2(x.size()); | |
171 | for (size_t i = 0; i < y.size(); ++i) | |
172 | { | |
173 | d2ydx2[i] = 6*x[i]; | |
174 | } | |
175 | ||
176 | auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)); | |
177 | ||
178 | for (Real t = 0; t <= 9; t += 0.0078125) | |
179 | { | |
180 | CHECK_ULP_CLOSE(t*t*t, qh(t), 10); | |
181 | CHECK_ULP_CLOSE(3*t*t, qh.prime(t), 15); | |
182 | CHECK_ULP_CLOSE(6*t, qh.double_prime(t), 20); | |
183 | } | |
184 | } | |
185 | ||
186 | template<typename Real> | |
187 | void test_quartic() | |
188 | { | |
189 | ||
190 | std::vector<Real> x{0,1,2,3, 4,5,6,7,8,9, 10, 11}; | |
191 | std::vector<Real> y(x.size()); | |
192 | for (size_t i = 0; i < y.size(); ++i) | |
193 | { | |
194 | y[i] = x[i]*x[i]*x[i]*x[i]; | |
195 | } | |
196 | ||
197 | std::vector<Real> dydx(x.size()); | |
198 | for (size_t i = 0; i < y.size(); ++i) | |
199 | { | |
200 | dydx[i] = 4*x[i]*x[i]*x[i]; | |
201 | } | |
202 | ||
203 | std::vector<Real> d2ydx2(x.size()); | |
204 | for (size_t i = 0; i < y.size(); ++i) | |
205 | { | |
206 | d2ydx2[i] = 12*x[i]*x[i]; | |
207 | } | |
208 | ||
209 | auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)); | |
210 | ||
211 | for (Real t = 1; t <= 11; t += 0.0078125) | |
212 | { | |
213 | CHECK_ULP_CLOSE(t*t*t*t, qh(t), 100); | |
214 | CHECK_ULP_CLOSE(4*t*t*t, qh.prime(t), 100); | |
215 | CHECK_ULP_CLOSE(12*t*t, qh.double_prime(t), 100); | |
216 | } | |
217 | } | |
218 | ||
219 | ||
220 | template<typename Real> | |
221 | void test_interpolation_condition() | |
222 | { | |
223 | for (size_t n = 4; n < 50; ++n) { | |
224 | std::vector<Real> x(n); | |
225 | std::vector<Real> y(n); | |
226 | std::vector<Real> dydx(n); | |
227 | std::vector<Real> d2ydx2(n); | |
228 | boost::random::mt19937 rd; | |
229 | boost::random::uniform_real_distribution<Real> dis(0,1); | |
230 | Real x0 = dis(rd); | |
231 | x[0] = x0; | |
232 | y[0] = dis(rd); | |
233 | for (size_t i = 1; i < n; ++i) { | |
234 | x[i] = x[i-1] + dis(rd); | |
235 | y[i] = dis(rd); | |
236 | dydx[i] = dis(rd); | |
237 | d2ydx2[i] = dis(rd); | |
238 | } | |
239 | ||
240 | auto x_copy = x; | |
241 | auto y_copy = y; | |
242 | auto dydx_copy = dydx; | |
243 | auto d2ydx2_copy = d2ydx2; | |
244 | auto s = quintic_hermite(std::move(x_copy), std::move(y_copy), std::move(dydx_copy), std::move(d2ydx2_copy)); | |
245 | //std::cout << "s = " << s << "\n"; | |
246 | for (size_t i = 0; i < x.size(); ++i) { | |
247 | CHECK_ULP_CLOSE(y[i], s(x[i]), 2); | |
248 | CHECK_ULP_CLOSE(dydx[i], s.prime(x[i]), 2); | |
249 | CHECK_ULP_CLOSE(d2ydx2[i], s.double_prime(x[i]), 2); | |
250 | } | |
251 | } | |
252 | } | |
253 | ||
254 | template<typename Real> | |
255 | void test_cardinal_constant() | |
256 | { | |
257 | ||
258 | std::vector<Real> y(25); | |
259 | std::vector<Real> dydx(y.size(), 0); | |
260 | std::vector<Real> d2ydx2(y.size(), 0); | |
261 | for (auto & t : y) { | |
262 | t = 7; | |
263 | } | |
264 | Real x0 = 4; | |
265 | Real dx = Real(1)/Real(8); | |
266 | ||
267 | auto qh = cardinal_quintic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), x0, dx); | |
268 | ||
269 | for (Real t = x0; t <= x0 + 24*dx; t += 0.25) | |
270 | { | |
271 | CHECK_ULP_CLOSE(Real(7), qh(t), 24); | |
272 | CHECK_ULP_CLOSE(Real(0), qh.prime(t), 24); | |
273 | CHECK_ULP_CLOSE(Real(0), qh.double_prime(t), 24); | |
274 | } | |
275 | ||
276 | std::vector<std::array<Real, 3>> data(25); | |
277 | for (size_t i = 0; i < data.size(); ++i) | |
278 | { | |
279 | data[i][0] = 7; | |
280 | data[i][1] = 0; | |
281 | data[i][2] = 0; | |
282 | } | |
283 | ||
284 | auto qh_aos = cardinal_quintic_hermite_aos(std::move(data), x0, dx); | |
285 | for (Real t = x0; t <= x0 + 24*dx; t += 0.25) | |
286 | { | |
287 | CHECK_ULP_CLOSE(Real(7), qh_aos(t), 24); | |
288 | CHECK_ULP_CLOSE(Real(0), qh_aos.prime(t), 24); | |
289 | CHECK_ULP_CLOSE(Real(0), qh_aos.double_prime(t), 24); | |
290 | } | |
291 | ||
292 | // Now check the boundaries: | |
293 | auto [tlo, thi] = qh.domain(); | |
294 | int samples = 5000; | |
295 | int i = 0; | |
296 | while (i++ < samples) | |
297 | { | |
298 | CHECK_ULP_CLOSE(Real(7), qh(tlo), 2); | |
299 | CHECK_ULP_CLOSE(Real(7), qh(thi), 2); | |
300 | CHECK_ULP_CLOSE(Real(7), qh_aos(tlo), 2); | |
301 | CHECK_ULP_CLOSE(Real(7), qh_aos(thi), 2); | |
302 | CHECK_ULP_CLOSE(Real(0), qh.prime(tlo), 2); | |
303 | CHECK_ULP_CLOSE(Real(0), qh.prime(thi), 2); | |
304 | CHECK_ULP_CLOSE(Real(0), qh_aos.prime(tlo), 2); | |
305 | CHECK_ULP_CLOSE(Real(0), qh_aos.prime(thi), 2); | |
306 | ||
1e59de90 | 307 | tlo = boost::math::nextafter(tlo, (std::numeric_limits<Real>::max)()); |
f67539c2 TL |
308 | thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest()); |
309 | } | |
310 | } | |
311 | ||
312 | ||
313 | template<typename Real> | |
314 | void test_cardinal_linear() | |
315 | { | |
316 | std::vector<Real> y{0,1,2,3,4,5,6,7,8,9}; | |
317 | Real x0 = 0; | |
318 | Real dx = 1; | |
319 | std::vector<Real> dydx(y.size(), 1); | |
320 | std::vector<Real> d2ydx2(y.size(), 0); | |
321 | ||
322 | auto qh = cardinal_quintic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), x0, dx); | |
323 | ||
324 | for (Real t = 0; t <= 9; t += 0.25) { | |
325 | CHECK_ULP_CLOSE(Real(t), qh(t), 2); | |
326 | CHECK_ULP_CLOSE(Real(1), qh.prime(t), 2); | |
327 | CHECK_ULP_CLOSE(Real(0), qh.double_prime(t), 2); | |
328 | } | |
329 | ||
330 | std::vector<std::array<Real, 3>> data(10); | |
331 | for (size_t i = 0; i < data.size(); ++i) { | |
332 | data[i][0] = i; | |
333 | data[i][1] = 1; | |
334 | data[i][2] = 0; | |
335 | } | |
336 | ||
337 | auto qh_aos = cardinal_quintic_hermite_aos(std::move(data), x0, dx); | |
338 | ||
339 | for (Real t = 0; t <= 9; t += 0.25) { | |
340 | CHECK_ULP_CLOSE(Real(t), qh_aos(t), 2); | |
341 | CHECK_ULP_CLOSE(Real(1), qh_aos.prime(t), 2); | |
342 | CHECK_ULP_CLOSE(Real(0), qh_aos.double_prime(t), 2); | |
343 | } | |
344 | ||
345 | // Now check the boundaries: | |
346 | auto [tlo, thi] = qh.domain(); | |
347 | int samples = 5000; | |
348 | int i = 0; | |
349 | while (i++ < samples) | |
350 | { | |
351 | CHECK_ULP_CLOSE(Real(tlo), qh(tlo), 2); | |
352 | CHECK_ULP_CLOSE(Real(thi), qh(thi), 2); | |
353 | CHECK_ULP_CLOSE(Real(tlo), qh_aos(tlo), 2); | |
354 | CHECK_ULP_CLOSE(Real(thi), qh_aos(thi), 2); | |
355 | CHECK_ULP_CLOSE(Real(1), qh.prime(tlo), 2); | |
356 | CHECK_ULP_CLOSE(Real(1), qh.prime(thi), 128); | |
357 | CHECK_ULP_CLOSE(Real(1), qh_aos.prime(tlo), 2); | |
358 | CHECK_ULP_CLOSE(Real(1), qh_aos.prime(thi), 128); | |
359 | ||
1e59de90 | 360 | tlo = boost::math::nextafter(tlo, (std::numeric_limits<Real>::max)()); |
f67539c2 TL |
361 | thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest()); |
362 | } | |
363 | } | |
364 | ||
365 | template<typename Real> | |
366 | void test_cardinal_quadratic() | |
367 | { | |
368 | Real x0 = 0; | |
369 | Real dx = 1; | |
370 | std::vector<Real> y(10); | |
371 | for (size_t i = 0; i < y.size(); ++i) | |
372 | { | |
373 | y[i] = i*i/Real(2); | |
374 | } | |
375 | ||
376 | std::vector<Real> dydx(y.size()); | |
377 | for (size_t i = 0; i < y.size(); ++i) { | |
378 | dydx[i] = i; | |
379 | } | |
380 | ||
381 | std::vector<Real> d2ydx2(y.size(), 1); | |
382 | ||
383 | auto qh = cardinal_quintic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), x0, dx); | |
384 | ||
385 | for (Real t = 0; t <= 9; t += 0.0078125) { | |
386 | Real computed = qh(t); | |
387 | CHECK_ULP_CLOSE(Real(t*t)/2, computed, 2); | |
388 | CHECK_ULP_CLOSE(t, qh.prime(t), 15); | |
389 | CHECK_ULP_CLOSE(Real(1), qh.double_prime(t), 32); | |
390 | } | |
391 | ||
392 | std::vector<std::array<Real, 3>> data(10); | |
393 | for (size_t i = 0; i < data.size(); ++i) { | |
394 | data[i][0] = i*i/Real(2); | |
395 | data[i][1] = i; | |
396 | data[i][2] = 1; | |
397 | } | |
398 | auto qh_aos = cardinal_quintic_hermite_aos(std::move(data), x0, dx); | |
399 | ||
400 | for (Real t = 0; t <= 9; t += 0.0078125) | |
401 | { | |
402 | Real computed = qh_aos(t); | |
403 | CHECK_ULP_CLOSE(Real(t*t)/2, computed, 2); | |
404 | CHECK_ULP_CLOSE(t, qh_aos.prime(t), 12); | |
405 | CHECK_ULP_CLOSE(Real(1), qh_aos.double_prime(t), 64); | |
406 | } | |
407 | ||
408 | // Now check the boundaries: | |
409 | auto [tlo, thi] = qh.domain(); | |
410 | int samples = 5000; | |
411 | int i = 0; | |
412 | while (i++ < samples) | |
413 | { | |
414 | CHECK_ULP_CLOSE(tlo*tlo/2, qh(tlo), 16); | |
415 | CHECK_ULP_CLOSE(thi*thi/2, qh(thi), 16); | |
416 | CHECK_ULP_CLOSE(tlo*tlo/2, qh_aos(tlo), 16); | |
417 | CHECK_ULP_CLOSE(thi*thi/2, qh_aos(thi), 16); | |
418 | CHECK_ULP_CLOSE(tlo, qh.prime(tlo), 16); | |
419 | CHECK_ULP_CLOSE(thi, qh.prime(thi), 64); | |
420 | CHECK_ULP_CLOSE(tlo, qh_aos.prime(tlo), 16); | |
421 | CHECK_ULP_CLOSE(thi, qh_aos.prime(thi), 64); | |
422 | ||
1e59de90 | 423 | tlo = boost::math::nextafter(tlo, (std::numeric_limits<Real>::max)()); |
f67539c2 TL |
424 | thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest()); |
425 | } | |
426 | } | |
427 | ||
428 | template<typename Real> | |
429 | void test_cardinal_cubic() | |
430 | { | |
431 | Real x0 = 0; | |
432 | Real dx = 1; | |
433 | std::vector<Real> y(10); | |
434 | for (size_t i = 0; i < y.size(); ++i) | |
435 | { | |
436 | y[i] = i*i*i; | |
437 | } | |
438 | ||
439 | std::vector<Real> dydx(y.size()); | |
440 | for (size_t i = 0; i < y.size(); ++i) { | |
441 | dydx[i] = 3*i*i; | |
442 | } | |
443 | ||
444 | std::vector<Real> d2ydx2(y.size()); | |
445 | for (size_t i = 0; i < y.size(); ++i) { | |
446 | d2ydx2[i] = 6*i; | |
447 | } | |
448 | ||
449 | auto qh = cardinal_quintic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), x0, dx); | |
450 | ||
451 | for (Real t = 0; t <= 9; t += 0.0078125) | |
452 | { | |
453 | Real computed = qh(t); | |
454 | CHECK_ULP_CLOSE(t*t*t, computed, 10); | |
455 | CHECK_ULP_CLOSE(3*t*t, qh.prime(t), 15); | |
456 | CHECK_ULP_CLOSE(6*t, qh.double_prime(t), 39); | |
457 | } | |
458 | ||
459 | std::vector<std::array<Real, 3>> data(10); | |
460 | for (size_t i = 0; i < data.size(); ++i) { | |
461 | data[i][0] = i*i*i; | |
462 | data[i][1] = 3*i*i; | |
463 | data[i][2] = 6*i; | |
464 | } | |
465 | ||
466 | auto qh_aos = cardinal_quintic_hermite_aos(std::move(data), x0, dx); | |
467 | for (Real t = 0; t <= 9; t += 0.0078125) | |
468 | { | |
469 | Real computed = qh_aos(t); | |
470 | CHECK_ULP_CLOSE(t*t*t, computed, 10); | |
471 | CHECK_ULP_CLOSE(3*t*t, qh_aos.prime(t), 15); | |
472 | CHECK_ULP_CLOSE(6*t, qh_aos.double_prime(t), 30); | |
473 | } | |
474 | } | |
475 | ||
476 | template<typename Real> | |
477 | void test_cardinal_quartic() | |
478 | { | |
479 | Real x0 = 0; | |
480 | Real dx = 1; | |
481 | std::vector<Real> y(7); | |
482 | for (size_t i = 0; i < y.size(); ++i) | |
483 | { | |
484 | y[i] = i*i*i*i; | |
485 | } | |
486 | ||
487 | std::vector<Real> dydx(y.size()); | |
488 | for (size_t i = 0; i < y.size(); ++i) { | |
489 | dydx[i] = 4*i*i*i; | |
490 | } | |
491 | ||
492 | std::vector<Real> d2ydx2(y.size()); | |
493 | for (size_t i = 0; i < y.size(); ++i) { | |
494 | d2ydx2[i] = 12*i*i; | |
495 | } | |
496 | ||
497 | auto qh = cardinal_quintic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), x0, dx); | |
498 | ||
499 | for (Real t = 0; t <= 6; t += 0.0078125) | |
500 | { | |
501 | CHECK_ULP_CLOSE(Real(t*t*t*t), qh(t), 250); | |
502 | CHECK_ULP_CLOSE(4*t*t*t, qh.prime(t), 250); | |
503 | CHECK_ULP_CLOSE(12*t*t, qh.double_prime(t), 250); | |
504 | } | |
505 | ||
506 | std::vector<std::array<Real, 3>> data(7); | |
507 | for (size_t i = 0; i < data.size(); ++i) { | |
508 | data[i][0] = i*i*i*i; | |
509 | data[i][1] = 4*i*i*i; | |
510 | data[i][2] = 12*i*i; | |
511 | } | |
512 | ||
513 | auto qh_aos = cardinal_quintic_hermite_aos(std::move(data), x0, dx); | |
514 | for (Real t = 0; t <= 6; t += 0.0078125) | |
515 | { | |
516 | Real computed = qh_aos(t); | |
517 | CHECK_ULP_CLOSE(t*t*t*t, computed, 10); | |
518 | CHECK_ULP_CLOSE(4*t*t*t, qh_aos.prime(t), 64); | |
519 | CHECK_ULP_CLOSE(12*t*t, qh_aos.double_prime(t), 128); | |
520 | } | |
521 | } | |
522 | ||
523 | ||
524 | int main() | |
525 | { | |
526 | test_constant<float>(); | |
527 | test_linear<float>(); | |
528 | test_quadratic<float>(); | |
529 | test_cubic<float>(); | |
530 | test_quartic<float>(); | |
531 | test_interpolation_condition<float>(); | |
532 | ||
533 | test_cardinal_constant<float>(); | |
534 | test_cardinal_linear<float>(); | |
535 | test_cardinal_quadratic<float>(); | |
536 | test_cardinal_cubic<float>(); | |
537 | test_cardinal_quartic<float>(); | |
538 | ||
539 | test_constant<double>(); | |
540 | test_linear<double>(); | |
541 | test_quadratic<double>(); | |
542 | test_cubic<double>(); | |
543 | test_quartic<double>(); | |
544 | test_interpolation_condition<double>(); | |
545 | ||
546 | test_cardinal_constant<double>(); | |
547 | test_cardinal_linear<double>(); | |
548 | test_cardinal_quadratic<double>(); | |
549 | test_cardinal_cubic<double>(); | |
550 | test_cardinal_quartic<double>(); | |
551 | ||
552 | test_constant<long double>(); | |
553 | test_linear<long double>(); | |
554 | test_quadratic<long double>(); | |
555 | test_cubic<long double>(); | |
556 | test_quartic<long double>(); | |
557 | test_interpolation_condition<long double>(); | |
558 | ||
559 | test_cardinal_constant<long double>(); | |
560 | test_cardinal_linear<long double>(); | |
561 | test_cardinal_quadratic<long double>(); | |
562 | test_cardinal_cubic<long double>(); | |
563 | test_cardinal_quartic<long double>(); | |
564 | ||
565 | #ifdef BOOST_HAS_FLOAT128 | |
566 | test_constant<float128>(); | |
567 | //test_linear<float128>(); | |
568 | test_quadratic<float128>(); | |
569 | test_cubic<float128>(); | |
570 | test_quartic<float128>(); | |
571 | test_interpolation_condition<float128>(); | |
572 | test_cardinal_constant<float128>(); | |
573 | test_cardinal_linear<float128>(); | |
574 | test_cardinal_quadratic<float128>(); | |
575 | test_cardinal_cubic<float128>(); | |
576 | test_cardinal_quartic<float128>(); | |
577 | #endif | |
578 | ||
579 | return boost::math::test::report_errors(); | |
580 | } |