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