]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
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())) { | |
f67539c2 | 172 | std::cerr << " Problem occurred at abscissa " << t << "\n"; |
92f5a8d4 TL |
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())) { | |
f67539c2 | 178 | std::cerr << " Problem occurred at abscissa " << t << "\n"; |
92f5a8d4 TL |
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())) { | |
f67539c2 | 184 | std::cerr << " Problem occurred at abscissa " << t << "\n"; |
92f5a8d4 TL |
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 | } |