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>
test_constant()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>
test_linear()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>
test_quadratic()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>
test_cubic()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>
test_quartic()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>
test_interpolation_condition()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>
test_cardinal_constant()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
307 tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
308 thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
309 }
310 }
311
312
313 template<typename Real>
test_cardinal_linear()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
360 tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
361 thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
362 }
363 }
364
365 template<typename Real>
test_cardinal_quadratic()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
423 tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
424 thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
425 }
426 }
427
428 template<typename Real>
test_cardinal_cubic()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>
test_cardinal_quartic()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
main()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 }
581