diff --git a/doc/interpolators/whittaker_shannon.qbk b/doc/interpolators/whittaker_shannon.qbk
index c025f644b..701651f53 100644
--- a/doc/interpolators/whittaker_shannon.qbk
+++ b/doc/interpolators/whittaker_shannon.qbk
@@ -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]
diff --git a/include/boost/math/interpolators/detail/whittaker_shannon_detail.hpp b/include/boost/math/interpolators/detail/whittaker_shannon_detail.hpp
index 2acfcb16c..3d4ace397 100644
--- a/include/boost/math/interpolators/detail/whittaker_shannon_detail.hpp
+++ b/include/boost/math/interpolators/detail/whittaker_shannon_detail.hpp
@@ -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;
     }
 
 
diff --git a/test/whittaker_shannon_test.cpp b/test/whittaker_shannon_test.cpp
index 3f61cc7c1..7bf24647b 100644
--- a/test/whittaker_shannon_test.cpp
+++ b/test/whittaker_shannon_test.cpp
@@ -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>();