123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233 |
- /*
- * Copyright Nick Thompson, 2019
- * Use, modification and distribution are subject to the
- * Boost Software License, Version 1.0. (See accompanying file
- * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
- */
- #include "math_unit_test.hpp"
- #include <vector>
- #include <random>
- #include <boost/math/constants/constants.hpp>
- #include <boost/math/interpolators/cardinal_trigonometric.hpp>
- #ifdef BOOST_HAS_FLOAT128
- #include <boost/multiprecision/float128.hpp>
- #endif
- using std::sin;
- using boost::math::constants::two_pi;
- using boost::math::interpolators::cardinal_trigonometric;
- template<class Real>
- void test_constant()
- {
- Real t0 = 0;
- Real h = 1;
- for(size_t n = 1; n < 20; ++n)
- {
- Real c = 8;
- std::vector<Real> v(n, c);
- auto ct = cardinal_trigonometric<decltype(v)>(v, t0, h);
- CHECK_ULP_CLOSE(c, ct(0.3), 3);
- CHECK_ULP_CLOSE(c*h*n, ct.integrate(), 3);
- CHECK_ULP_CLOSE(c*c*h*n, ct.squared_l2(), 3);
- CHECK_MOLLIFIED_CLOSE(Real(0), ct.prime(0.8), 25*std::numeric_limits<Real>::epsilon());
- CHECK_MOLLIFIED_CLOSE(Real(0), ct.double_prime(0.8), 25*std::numeric_limits<Real>::epsilon());
- }
- }
- template<class Real>
- void test_interpolation_condition()
- {
- std::mt19937 gen(1234);
- std::uniform_real_distribution<Real> dis(1, 10);
- for(size_t n = 1; n < 20; ++n) {
- Real t0 = dis(gen);
- Real h = dis(gen);
- std::vector<Real> v(n);
- for (size_t i = 0; i < n; ++i) {
- v[i] = dis(gen);
- }
- auto ct = cardinal_trigonometric<decltype(v)>(v, t0, h);
- for (size_t i = 0; i < n; ++i) {
- Real arg = t0 + i*h;
- Real expected = v[i];
- Real computed = ct(arg);
- if(!CHECK_ULP_CLOSE(expected, computed, 5*n))
- {
- std::cerr << " Samples: " << n << "\n";
- }
- }
- }
- }
- #ifdef BOOST_HAS_FLOAT128
- void test_constant_q()
- {
- __float128 t0 = 0;
- __float128 h = 1;
- for(size_t n = 1; n < 20; ++n)
- {
- __float128 c = 8;
- std::vector<__float128> v(n, c);
- auto ct = cardinal_trigonometric<decltype(v)>(v, t0, h);
- CHECK_ULP_CLOSE(boost::multiprecision::float128(c), boost::multiprecision::float128(ct(0.3)), 3);
- CHECK_ULP_CLOSE(boost::multiprecision::float128(c*h*n), boost::multiprecision::float128(ct.integrate()), 3);
- }
- }
- #endif
- template<class Real>
- void test_sampled_sine()
- {
- using std::sin;
- using std::cos;
- for (unsigned n = 15; n < 50; ++n)
- {
- Real t0 = 0;
- Real T = 1;
- Real h = T/n;
- std::vector<Real> v(n);
- auto s = [&](Real t) { return sin(two_pi<Real>()*(t-t0)/T);};
- auto s_prime = [&](Real t) { return two_pi<Real>()*cos(two_pi<Real>()*(t-t0)/T)/T;};
- auto s_double_prime = [&](Real t) { return -two_pi<Real>()*two_pi<Real>()*sin(two_pi<Real>()*(t-t0)/T)/(T*T);};
- for(size_t j = 0; j < v.size(); ++j)
- {
- Real t = t0 + j*h;
- v[j] = s(t);
- }
- auto ct = cardinal_trigonometric<decltype(v)>(v, t0, h);
- CHECK_ULP_CLOSE(T, ct.period(), 3);
- std::mt19937 gen(1234);
- std::uniform_real_distribution<Real> dist(0, 500);
- unsigned j = 0;
- while (j++ < 50) {
- Real arg = dist(gen);
- Real expected = s(arg);
- Real computed = ct(arg);
- CHECK_MOLLIFIED_CLOSE(expected, computed, std::numeric_limits<Real>::epsilon()*4000);
- expected = s_prime(arg);
- computed = ct.prime(arg);
- CHECK_MOLLIFIED_CLOSE(expected, computed, 18000*std::numeric_limits<Real>::epsilon());
- expected = s_double_prime(arg);
- computed = ct.double_prime(arg);
- CHECK_MOLLIFIED_CLOSE(expected, computed, 100000*std::numeric_limits<Real>::epsilon());
- }
- CHECK_MOLLIFIED_CLOSE(Real(0), ct.integrate(), std::numeric_limits<Real>::epsilon());
- }
- }
- template<class Real>
- void test_bump()
- {
- using std::exp;
- using std::abs;
- using std::sqrt;
- using std::pow;
- auto bump = [](Real x)->Real { if (abs(x) >= 1) { return Real(0); } return exp(-Real(1)/(Real(1)-x*x)); };
- auto bump_prime = [](Real x)->Real {
- if (abs(x) >= 1) { return Real(0); }
- return -2*x*exp(-Real(1)/(Real(1)-x*x))/pow(1-x*x,2);
- };
- auto bump_double_prime = [](Real x)->Real {
- if (abs(x) >= 1) { return Real(0); }
- return (6*pow(x,4)-2)*exp(-Real(1)/(Real(1)-x*x))/pow(1-x*x,4);
- };
- Real t0 = -1;
- size_t n = 4096;
- Real h = Real(2)/Real(n);
- std::vector<Real> v(n);
- for(size_t i = 0; i < n; ++i)
- {
- Real t = t0 + i*h;
- v[i] = bump(t);
- }
- auto ct = cardinal_trigonometric<decltype(v)>(v, t0, h);
- std::mt19937 gen(323723);
- std::uniform_real_distribution<long double> dis(-0.9, 0.9);
- size_t i = 0;
- while (i++ < 1000)
- {
- Real t = static_cast<Real>(dis(gen));
- Real expected = bump(t);
- Real computed = ct(t);
- if(!CHECK_MOLLIFIED_CLOSE(expected, computed, 2*std::numeric_limits<Real>::epsilon())) {
- std::cerr << " Problem occured at abscissa " << t << "\n";
- }
- expected = bump_prime(t);
- computed = ct.prime(t);
- if(!CHECK_MOLLIFIED_CLOSE(expected, computed, 4000*std::numeric_limits<Real>::epsilon())) {
- std::cerr << " Problem occured at abscissa " << t << "\n";
- }
- expected = bump_double_prime(t);
- computed = ct.double_prime(t);
- if(!CHECK_MOLLIFIED_CLOSE(expected, computed, 4000*4000*std::numeric_limits<Real>::epsilon())) {
- std::cerr << " Problem occured at abscissa " << t << "\n";
- }
- }
- // Wolfram Alpha:
- // NIntegrate[Exp[-1/(1-x*x)],{x,-1,1}]
- CHECK_ULP_CLOSE(Real(0.443993816168079437823L), ct.integrate(), 3);
- // NIntegrate[Exp[-2/(1-x*x)],{x,-1,1}]
- CHECK_ULP_CLOSE(Real(0.1330861208449942715569473279553285713625791551628130055345002588895389L), ct.squared_l2(), 1);
- }
- int main()
- {
- #ifdef TEST1
- test_constant<float>();
- test_sampled_sine<float>();
- test_bump<float>();
- test_interpolation_condition<float>();
- #endif
- #ifdef TEST2
- test_constant<double>();
- test_sampled_sine<double>();
- test_bump<double>();
- test_interpolation_condition<double>();
- #endif
- #ifdef TEST3
- test_constant<long double>();
- test_sampled_sine<long double>();
- test_bump<long double>();
- test_interpolation_condition<long double>();
- #endif
- #ifdef TEST4
- #ifdef BOOST_HAS_FLOAT128
- test_constant_q();
- #endif
- #endif
- return boost::math::test::report_errors();
- }
|