boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Whittaker-Shannon: Update docs, implement derivatives.

Browse Source
This commit is contained in:
Nick Thompson
2019-06-24 08:28:14 -04:00
parent b2c2f0e644
commit 260a3af015
3 changed files with 17 additions and 27 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
doc/interpolators/whittaker_shannon.qbk
View File

@@ -24,6 +24,8 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
whittaker_shannon(RandomAccessContainer&& v, Real left_endpoint, Real step_size);
Real operator()(Real x) const;
Real prime(Real x) const;
};
}}} // namespaces
@@ -55,10 +57,13 @@ Here is an example of interpolating a smooth "bump function":
double y = ws(0.3);
The derivative of the interpolant can also be evaluated, but the accuracy is not as high:
double yp = ws.prime(0.3);
[heading Complexity and Performance]
The call to the constructor requires [bigo](1) operations, simply moving data into the class.
Each call the the interpolant is [bigo](/n/), where /n/ is the number of points to interpolate.
The compiler flag `-ffast-math` allows vectorization of the call to the interpolant (at least on clang) and is highly recommended.
[endsect] [/section:whittaker_shannon]

9
include/boost/math/interpolators/detail/whittaker_shannon_detail.hpp
View File

@@ -64,7 +64,6 @@ public:
Real x = (t - m_t0)/m_h;
if (ceil(x) == x) {
//std::cout << "TAKING INTEGER ARGUMENT BRANCH: x = " << x << "\n";
Real s = 0;
long j = static_cast<long>(x);
long n = m_y.size();
@@ -86,19 +85,17 @@ public:
Real z = x;
auto it = m_y.begin();
Real cospix = boost::math::cos_pi(x);
Real sinpix = boost::math::sin_pi(x);
Real sinpix_div_pi = boost::math::sin_pi(x)/pi<Real>();
Real s = 0;
// For some reason, neither clang nor g++ will cache the address of m_y.end() in a register.
// Hence make a copy of it:
auto end = m_y.end();
while(it != end)
{
s += (*it++)*(pi<Real>()*z*cospix - sinpix)/(z*z);
s += (*it++)*(z*cospix - sinpix_div_pi)/(z*z);
z -= 1;
}
return s/(pi<Real>()*m_h);
return s/m_h;
}

28
test/whittaker_shannon_test.cpp
View File

@@ -37,7 +37,6 @@ void test_trivial()
if(!CHECK_MOLLIFIED_CLOSE(expected, ws.prime(h), 10*std::numeric_limits<Real>::epsilon())) {
std::cerr << " Problem occured at abscissa " << 0 << "\n";
}
}
template<class Real>
@@ -48,12 +47,11 @@ void test_knots()
size_t n = 512;
std::vector<Real> v(n);
std::mt19937 gen(323723);
std::uniform_real_distribution<long double> dis(1.0, 2.0);
std::uniform_real_distribution<Real> dis(1.0, 2.0);
for(size_t i = 0; i < n; ++i) {
v[i] = static_cast<Real>(dis(gen));
}
auto ws = whittaker_shannon<decltype(v)>(std::move(v), t0, h);
size_t i = 0;
@@ -61,14 +59,7 @@ void test_knots()
Real t = t0 + i*h;
Real expected = ws[i];
Real computed = ws(t);
using std::isnan;
if(isnan(computed)) {
throw std::domain_error("This is bad m'kay?");
}
if(isnan(expected)) {
throw std::domain_error("This is bad m'kay?");
}
CHECK_ULP_CLOSE(ws[i], ws(t), 16);
CHECK_ULP_CLOSE(expected, computed, 16);
++i;
}
}
@@ -78,6 +69,7 @@ void test_bump()
{
using std::exp;
using std::abs;
using std::sqrt;
auto bump = [](Real x) { if (abs(x) >= 1) { return Real(0); } return exp(-Real(1)/(Real(1)-x*x)); };
auto bump_prime = [&bump](Real x) { Real z = 1-x*x; return -2*x*bump(x)/(z*z); };
@@ -107,14 +99,14 @@ void test_bump()
Real expected_prime = bump_prime(t);
Real computed_prime = ws.prime(t);
if(!CHECK_MOLLIFIED_CLOSE(expected_prime, computed_prime, 10*std::numeric_limits<Real>::epsilon())) {
if(!CHECK_MOLLIFIED_CLOSE(expected_prime, computed_prime, 1000*std::numeric_limits<Real>::epsilon())) {
std::cerr << " Problem occured at abscissa " << t << "\n";
}
}
std::mt19937 gen(323723);
std::uniform_real_distribution<long double> dis(-0.95, 0.95);
std::uniform_real_distribution<long double> dis(-0.85, 0.85);
size_t i = 0;
while (i++ < 1000)
@@ -128,7 +120,7 @@ void test_bump()
Real expected_prime = bump_prime(t);
Real computed_prime = ws.prime(t);
if(!CHECK_MOLLIFIED_CLOSE(expected_prime, computed_prime, 10*std::numeric_limits<Real>::epsilon())) {
if(!CHECK_MOLLIFIED_CLOSE(expected_prime, computed_prime, sqrt(std::numeric_limits<Real>::epsilon()))) {
std::cerr << " Problem occured at abscissa " << t << "\n";
}
}
@@ -137,16 +129,12 @@ void test_bump()
int main()
{
/* test_knots<float>();
test_knots<float>();
test_knots<double>();
test_knots<long double>();
#ifdef BOOST_HAS_FLOAT128
test_knots<float128>();
#endif*/
//test_bump<float>();
test_bump<double>();
//test_bump<long double>();
test_bump<long double>();
test_trivial<float>();
test_trivial<double>();