libs/math/test/cardinal_trigonometric_test.cpp

   1 /*
   2  * Copyright Nick Thompson, 2019
   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 <vector>
  10 #include <random>
  11 #include <boost/math/constants/constants.hpp>
  12 #include <boost/math/interpolators/cardinal_trigonometric.hpp>
  13 #ifdef BOOST_HAS_FLOAT128
  14 #include <boost/multiprecision/float128.hpp>
  15 #endif
  16
  17
  18 using std::sin;
  19 using boost::math::constants::two_pi;
  20 using boost::math::interpolators::cardinal_trigonometric;
  21
  22 template<class Real>
  23 void test_constant()
  24 {
  25     Real t0 = 0;
  26     Real h = 1;
  27     for(size_t n = 1; n < 20; ++n)
  28     {
  29       Real c = 8;
  30       std::vector<Real> v(n, c);
  31       auto ct = cardinal_trigonometric<decltype(v)>(v, t0, h);
  32       CHECK_ULP_CLOSE(c, ct(0.3), 3);
  33       CHECK_ULP_CLOSE(c*h*n, ct.integrate(), 3);
  34       CHECK_ULP_CLOSE(c*c*h*n, ct.squared_l2(), 3);
  35       CHECK_MOLLIFIED_CLOSE(Real(0), ct.prime(0.8), 25*std::numeric_limits<Real>::epsilon());
  36       CHECK_MOLLIFIED_CLOSE(Real(0), ct.double_prime(0.8), 25*std::numeric_limits<Real>::epsilon());
  37     }
  38 }
  39
  40 template<class Real>
  41 void test_interpolation_condition()
  42 {
  43   std::mt19937 gen(1234);
  44   std::uniform_real_distribution<Real> dis(1, 10);
  45
  46   for(size_t n = 1; n < 20; ++n) {
  47     Real t0 = dis(gen);
  48     Real h = dis(gen);
  49     std::vector<Real> v(n);
  50     for (size_t i = 0; i < n; ++i) {
  51       v[i] = dis(gen);
  52     }
  53     auto ct = cardinal_trigonometric<decltype(v)>(v, t0, h);
  54     for (size_t i = 0; i < n; ++i) {
  55       Real arg = t0 + i*h;
  56       Real expected = v[i];
  57       Real computed = ct(arg);
  58       if(!CHECK_ULP_CLOSE(expected, computed, 5*n))
  59       {
  60         std::cerr << "  Samples: " << n << "\n";
  61       }
  62     }
  63   }
  64
  65 }
  66
  67
  68 #ifdef BOOST_HAS_FLOAT128
  69 void test_constant_q()
  70 {
  71     __float128 t0 = 0;
  72     __float128 h = 1;
  73     for(size_t n = 1; n < 20; ++n)
  74     {
  75       __float128 c = 8;
  76       std::vector<__float128> v(n, c);
  77       auto ct = cardinal_trigonometric<decltype(v)>(v, t0, h);
  78       CHECK_ULP_CLOSE(boost::multiprecision::float128(c), boost::multiprecision::float128(ct(0.3)), 3);
  79       CHECK_ULP_CLOSE(boost::multiprecision::float128(c*h*n), boost::multiprecision::float128(ct.integrate()), 3);
  80     }
  81 }
  82 #endif
  83
  84
  85 template<class Real>
  86 void test_sampled_sine()
  87 {
  88     using std::sin;
  89     using std::cos;
  90     for (unsigned n = 15; n < 50; ++n)
  91     {
  92       Real t0 = 0;
  93       Real T = 1;
  94       Real h = T/n;
  95       std::vector<Real> v(n);
  96       auto s = [&](Real t) { return sin(two_pi<Real>()*(t-t0)/T);};
  97       auto s_prime = [&](Real t) { return two_pi<Real>()*cos(two_pi<Real>()*(t-t0)/T)/T;};
  98       auto s_double_prime = [&](Real t) { return -two_pi<Real>()*two_pi<Real>()*sin(two_pi<Real>()*(t-t0)/T)/(T*T);};
  99       for(size_t j = 0; j < v.size(); ++j)
 100       {
 101           Real t = t0 + j*h;
 102           v[j] = s(t);
 103       }
 104       auto ct = cardinal_trigonometric<decltype(v)>(v, t0, h);
 105       CHECK_ULP_CLOSE(T, ct.period(), 3);
 106       std::mt19937 gen(1234);
 107       std::uniform_real_distribution<Real> dist(0, 500);
 108
 109       unsigned j = 0;
 110       while (j++ < 50) {
 111         Real arg = dist(gen);
 112         Real expected = s(arg);
 113         Real computed = ct(arg);
 114         CHECK_MOLLIFIED_CLOSE(expected, computed, std::numeric_limits<Real>::epsilon()*4000);
 115
 116         expected = s_prime(arg);
 117         computed = ct.prime(arg);
 118         CHECK_MOLLIFIED_CLOSE(expected, computed, 18000*std::numeric_limits<Real>::epsilon());
 119
 120         expected = s_double_prime(arg);
 121         computed = ct.double_prime(arg);
 122         CHECK_MOLLIFIED_CLOSE(expected, computed, 100000*std::numeric_limits<Real>::epsilon());
 123
 124       }
 125       CHECK_MOLLIFIED_CLOSE(Real(0), ct.integrate(), std::numeric_limits<Real>::epsilon());
 126     }
 127 }
 128
 129 template<class Real>
 130 void test_bump()
 131 {
 132   using std::exp;
 133   using std::abs;
 134   using std::sqrt;
 135   using std::pow;
 136   auto bump = [](Real x)->Real { if (abs(x) >= 1) { return Real(0); } return exp(-Real(1)/(Real(1)-x*x)); };
 137   auto bump_prime = [](Real x)->Real {
 138       if (abs(x) >= 1) { return Real(0); }
 139
 140       return -2*x*exp(-Real(1)/(Real(1)-x*x))/pow(1-x*x,2);
 141   };
 142
 143   auto bump_double_prime = [](Real x)->Real {
 144       if (abs(x) >= 1) { return Real(0); }
 145
 146       return (6*pow(x,4)-2)*exp(-Real(1)/(Real(1)-x*x))/pow(1-x*x,4);
 147   };
 148
 149
 150   Real t0 = -1;
 151   size_t n = 4096;
 152   Real h = Real(2)/Real(n);
 153
 154   std::vector<Real> v(n);
 155   for(size_t i = 0; i < n; ++i)
 156   {
 157       Real t = t0 + i*h;
 158       v[i] = bump(t);
 159   }
 160
 161   auto ct = cardinal_trigonometric<decltype(v)>(v, t0, h);
 162   std::mt19937 gen(323723);
 163   std::uniform_real_distribution<long double> dis(-0.9, 0.9);
 164
 165   size_t i = 0;
 166   while (i++ < 1000)
 167   {
 168       Real t = static_cast<Real>(dis(gen));
 169       Real expected = bump(t);
 170       Real computed = ct(t);
 171       if(!CHECK_MOLLIFIED_CLOSE(expected, computed, 2*std::numeric_limits<Real>::epsilon())) {
 172           std::cerr << "  Problem occured at abscissa " << t << "\n";
 173       }
 174
 175       expected = bump_prime(t);
 176       computed = ct.prime(t);
 177       if(!CHECK_MOLLIFIED_CLOSE(expected, computed, 4000*std::numeric_limits<Real>::epsilon())) {
 178           std::cerr << "  Problem occured at abscissa " << t << "\n";
 179       }
 180
 181       expected = bump_double_prime(t);
 182       computed = ct.double_prime(t);
 183       if(!CHECK_MOLLIFIED_CLOSE(expected, computed, 4000*4000*std::numeric_limits<Real>::epsilon())) {
 184           std::cerr << "  Problem occured at abscissa " << t << "\n";
 185       }
 186
 187
 188   }
 189
 190   // Wolfram Alpha:
 191   // NIntegrate[Exp[-1/(1-x*x)],{x,-1,1}]
 192   CHECK_ULP_CLOSE(Real(0.443993816168079437823L), ct.integrate(), 3);
 193
 194   // NIntegrate[Exp[-2/(1-x*x)],{x,-1,1}]
 195   CHECK_ULP_CLOSE(Real(0.1330861208449942715569473279553285713625791551628130055345002588895389L), ct.squared_l2(), 1);
 196
 197
 198 }
 199
 200
 201 int main()
 202 {
 203
 204 #ifdef TEST1
 205     test_constant<float>();
 206     test_sampled_sine<float>();
 207     test_bump<float>();
 208     test_interpolation_condition<float>();
 209 #endif
 210
 211
 212 #ifdef TEST2
 213     test_constant<double>();
 214     test_sampled_sine<double>();
 215     test_bump<double>();
 216     test_interpolation_condition<double>();
 217 #endif
 218
 219 #ifdef TEST3
 220     test_constant<long double>();
 221     test_sampled_sine<long double>();
 222     test_bump<long double>();
 223     test_interpolation_condition<long double>();
 224 #endif
 225
 226 #ifdef TEST4
 227 #ifdef BOOST_HAS_FLOAT128
 228 test_constant_q();
 229 #endif
 230 #endif
 231
 232     return boost::math::test::report_errors();
 233 }