summaryrefslogtreecommitdiff
path: root/boost/math/interpolators
diff options
context:
space:
mode:
Diffstat (limited to 'boost/math/interpolators')
-rw-r--r--boost/math/interpolators/barycentric_rational.hpp70
-rw-r--r--boost/math/interpolators/cubic_b_spline.hpp78
-rw-r--r--boost/math/interpolators/detail/barycentric_rational_detail.hpp151
-rw-r--r--boost/math/interpolators/detail/cubic_b_spline_detail.hpp287
4 files changed, 586 insertions, 0 deletions
diff --git a/boost/math/interpolators/barycentric_rational.hpp b/boost/math/interpolators/barycentric_rational.hpp
new file mode 100644
index 0000000000..79bab9042d
--- /dev/null
+++ b/boost/math/interpolators/barycentric_rational.hpp
@@ -0,0 +1,70 @@
+/*
+ * Copyright Nick Thompson, 2017
+ * 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)
+ *
+ * Given N samples (t_i, y_i) which are irregularly spaced, this routine constructs an
+ * interpolant s which is constructed in O(N) time, occupies O(N) space, and can be evaluated in O(N) time.
+ * The interpolation is stable, unless one point is incredibly close to another, and the next point is incredibly far.
+ * The measure of this stability is the "local mesh ratio", which can be queried from the routine.
+ * Pictorially, the following t_i spacing is bad (has a high local mesh ratio)
+ * || | | | |
+ * and this t_i spacing is good (has a low local mesh ratio)
+ * | | | | | | | | | |
+ *
+ *
+ * If f is C^{d+2}, then the interpolant is O(h^(d+1)) accurate, where d is the interpolation order.
+ * A disadvantage of this interpolant is that it does not reproduce rational functions; for example, 1/(1+x^2) is not interpolated exactly.
+ *
+ * References:
+ * Floater, Michael S., and Kai Hormann. "Barycentric rational interpolation with no poles and high rates of approximation." Numerische Mathematik 107.2 (2007): 315-331.
+ * Press, William H., et al. "Numerical recipes third edition: the art of scientific computing." Cambridge University Press 32 (2007): 10013-2473.
+ */
+
+#ifndef BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_HPP
+#define BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_HPP
+
+#include <memory>
+#include <boost/math/interpolators/detail/barycentric_rational_detail.hpp>
+
+namespace boost{ namespace math{
+
+template<class Real>
+class barycentric_rational
+{
+public:
+ barycentric_rational(const Real* const x, const Real* const y, size_t n, size_t approximation_order = 3);
+
+ template <class InputIterator1, class InputIterator2>
+ barycentric_rational(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order = 3, typename boost::disable_if_c<boost::is_integral<InputIterator2>::value>::type* = 0);
+
+ Real operator()(Real x) const;
+
+private:
+ std::shared_ptr<detail::barycentric_rational_imp<Real>> m_imp;
+};
+
+template <class Real>
+barycentric_rational<Real>::barycentric_rational(const Real* const x, const Real* const y, size_t n, size_t approximation_order):
+ m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(x, x + n, y, approximation_order))
+{
+ return;
+}
+
+template <class Real>
+template <class InputIterator1, class InputIterator2>
+barycentric_rational<Real>::barycentric_rational(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order, typename boost::disable_if_c<boost::is_integral<InputIterator2>::value>::type*)
+ : m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(start_x, end_x, start_y, approximation_order))
+{
+}
+
+template<class Real>
+Real barycentric_rational<Real>::operator()(Real x) const
+{
+ return m_imp->operator()(x);
+}
+
+
+}}
+#endif
diff --git a/boost/math/interpolators/cubic_b_spline.hpp b/boost/math/interpolators/cubic_b_spline.hpp
new file mode 100644
index 0000000000..73ac1d0137
--- /dev/null
+++ b/boost/math/interpolators/cubic_b_spline.hpp
@@ -0,0 +1,78 @@
+// Copyright Nick Thompson, 2017
+// 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)
+
+// This implements the compactly supported cubic b spline algorithm described in
+// Kress, Rainer. "Numerical analysis, volume 181 of Graduate Texts in Mathematics." (1998).
+// Splines of compact support are faster to evaluate and are better conditioned than classical cubic splines.
+
+// Let f be the function we are trying to interpolate, and s be the interpolating spline.
+// The routine constructs the interpolant in O(N) time, and evaluating s at a point takes constant time.
+// The order of accuracy depends on the regularity of the f, however, assuming f is
+// four-times continuously differentiable, the error is of O(h^4).
+// In addition, we can differentiate the spline and obtain a good interpolant for f'.
+// The main restriction of this method is that the samples of f must be evenly spaced.
+// Look for barycentric rational interpolation for non-evenly sampled data.
+// Properties:
+// - s(x_j) = f(x_j)
+// - All cubic polynomials interpolated exactly
+
+#ifndef BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP
+#define BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP
+
+#include <boost/math/interpolators/detail/cubic_b_spline_detail.hpp>
+
+namespace boost{ namespace math{
+
+template <class Real>
+class cubic_b_spline
+{
+public:
+ // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
+ // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
+ template <class BidiIterator>
+ cubic_b_spline(const BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+ Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+ Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
+ cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
+ Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+ Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
+
+ cubic_b_spline() = default;
+ Real operator()(Real x) const;
+
+ Real prime(Real x) const;
+
+private:
+ std::shared_ptr<detail::cubic_b_spline_imp<Real>> m_imp;
+};
+
+template<class Real>
+cubic_b_spline<Real>::cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
+ Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cubic_b_spline_imp<Real>>(f, f + length, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
+{
+}
+
+template <class Real>
+template <class BidiIterator>
+cubic_b_spline<Real>::cubic_b_spline(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+ Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cubic_b_spline_imp<Real>>(f, end_p, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
+{
+}
+
+template<class Real>
+Real cubic_b_spline<Real>::operator()(Real x) const
+{
+ return m_imp->operator()(x);
+}
+
+template<class Real>
+Real cubic_b_spline<Real>::prime(Real x) const
+{
+ return m_imp->prime(x);
+}
+
+}}
+#endif
diff --git a/boost/math/interpolators/detail/barycentric_rational_detail.hpp b/boost/math/interpolators/detail/barycentric_rational_detail.hpp
new file mode 100644
index 0000000000..b853901d89
--- /dev/null
+++ b/boost/math/interpolators/detail/barycentric_rational_detail.hpp
@@ -0,0 +1,151 @@
+/*
+ * Copyright Nick Thompson, 2017
+ * 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)
+ */
+
+#ifndef BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP
+#define BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP
+
+#include <vector>
+#include <boost/lexical_cast.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/core/demangle.hpp>
+
+namespace boost{ namespace math{ namespace detail{
+
+template<class Real>
+class barycentric_rational_imp
+{
+public:
+ template <class InputIterator1, class InputIterator2>
+ barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order = 3);
+
+ Real operator()(Real x) const;
+
+ // The barycentric weights are not really that interesting; except to the unit tests!
+ Real weight(size_t i) const { return m_w[i]; }
+
+private:
+ // Technically, we do not need to copy over y to m_y, or x to m_x.
+ // We could simply store a pointer. However, in doing so,
+ // we'd need to make sure the user managed the lifetime of m_y,
+ // and didn't alter its data. Since we are unlikely to run out of
+ // memory for a linearly scaling algorithm, it seems best just to make a copy.
+ std::vector<Real> m_y;
+ std::vector<Real> m_x;
+ std::vector<Real> m_w;
+};
+
+template <class Real>
+template <class InputIterator1, class InputIterator2>
+barycentric_rational_imp<Real>::barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order)
+{
+ using std::abs;
+
+ std::ptrdiff_t n = std::distance(start_x, end_x);
+
+ if (approximation_order >= (std::size_t)n)
+ {
+ throw std::domain_error("Approximation order must be < data length.");
+ }
+
+ // Big sad memcpy to make sure the object is easy to use.
+ m_x.resize(n);
+ m_y.resize(n);
+ for(unsigned i = 0; start_x != end_x; ++start_x, ++start_y, ++i)
+ {
+ // But if we're going to do a memcpy, we can do some error checking which is inexpensive relative to the copy:
+ if(boost::math::isnan(*start_x))
+ {
+ std::string msg = std::string("x[") + boost::lexical_cast<std::string>(i) + "] is a NAN";
+ throw std::domain_error(msg);
+ }
+
+ if(boost::math::isnan(*start_y))
+ {
+ std::string msg = std::string("y[") + boost::lexical_cast<std::string>(i) + "] is a NAN";
+ throw std::domain_error(msg);
+ }
+
+ m_x[i] = *start_x;
+ m_y[i] = *start_y;
+ }
+
+ m_w.resize(n, 0);
+ for(int64_t k = 0; k < n; ++k)
+ {
+ int64_t i_min = (std::max)(k - (int64_t) approximation_order, (int64_t) 0);
+ int64_t i_max = k;
+ if (k >= n - (std::ptrdiff_t)approximation_order)
+ {
+ i_max = n - approximation_order - 1;
+ }
+
+ for(int64_t i = i_min; i <= i_max; ++i)
+ {
+ Real inv_product = 1;
+ int64_t j_max = (std::min)(static_cast<int64_t>(i + approximation_order), static_cast<int64_t>(n - 1));
+ for(int64_t j = i; j <= j_max; ++j)
+ {
+ if (j == k)
+ {
+ continue;
+ }
+
+ Real diff = m_x[k] - m_x[j];
+ if (abs(diff) < std::numeric_limits<Real>::epsilon())
+ {
+ std::string msg = std::string("Spacing between x[")
+ + boost::lexical_cast<std::string>(k) + std::string("] and x[")
+ + boost::lexical_cast<std::string>(i) + std::string("] is ")
+ + boost::lexical_cast<std::string>(diff) + std::string(", which is smaller than the epsilon of ")
+ + boost::core::demangle(typeid(Real).name());
+ throw std::logic_error(msg);
+ }
+ inv_product *= diff;
+ }
+ if (i % 2 == 0)
+ {
+ m_w[k] += 1/inv_product;
+ }
+ else
+ {
+ m_w[k] -= 1/inv_product;
+ }
+ }
+ }
+}
+
+template<class Real>
+Real barycentric_rational_imp<Real>::operator()(Real x) const
+{
+ Real numerator = 0;
+ Real denominator = 0;
+ for(size_t i = 0; i < m_x.size(); ++i)
+ {
+ // Presumably we should see if the accuracy is improved by using ULP distance of say, 5 here, instead of testing for floating point equality.
+ // However, it has been shown that if x approx x_i, but x != x_i, then inaccuracy in the numerator cancels the inaccuracy in the denominator,
+ // and the result is fairly accurate. See: http://epubs.siam.org/doi/pdf/10.1137/S0036144502417715
+ if (x == m_x[i])
+ {
+ return m_y[i];
+ }
+ Real t = m_w[i]/(x - m_x[i]);
+ numerator += t*m_y[i];
+ denominator += t;
+ }
+ return numerator/denominator;
+}
+
+/*
+ * A formula for computing the derivative of the barycentric representation is given in
+ * "Some New Aspects of Rational Interpolation", by Claus Schneider and Wilhelm Werner,
+ * Mathematics of Computation, v47, number 175, 1986.
+ * http://www.ams.org/journals/mcom/1986-47-175/S0025-5718-1986-0842136-8/S0025-5718-1986-0842136-8.pdf
+ * However, this requires a lot of machinery which is not built into the library at present.
+ * So we wait until there is a requirement to interpolate the derivative.
+ */
+}}}
+#endif
diff --git a/boost/math/interpolators/detail/cubic_b_spline_detail.hpp b/boost/math/interpolators/detail/cubic_b_spline_detail.hpp
new file mode 100644
index 0000000000..f7b2d6cd29
--- /dev/null
+++ b/boost/math/interpolators/detail/cubic_b_spline_detail.hpp
@@ -0,0 +1,287 @@
+// Copyright Nick Thompson, 2017
+// 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)
+
+#ifndef CUBIC_B_SPLINE_DETAIL_HPP
+#define CUBIC_B_SPLINE_DETAIL_HPP
+
+#include <limits>
+#include <cmath>
+#include <vector>
+#include <memory>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+
+namespace boost{ namespace math{ namespace detail{
+
+
+template <class Real>
+class cubic_b_spline_imp
+{
+public:
+ // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
+ // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
+ template <class BidiIterator>
+ cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+ Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+ Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
+
+ Real operator()(Real x) const;
+
+ Real prime(Real x) const;
+
+private:
+ std::vector<Real> m_beta;
+ Real m_h_inv;
+ Real m_a;
+ Real m_avg;
+};
+
+
+
+template <class Real>
+Real b3_spline(Real x)
+{
+ using std::abs;
+ Real absx = abs(x);
+ if (absx < 1)
+ {
+ Real y = 2 - absx;
+ Real z = 1 - absx;
+ return boost::math::constants::sixth<Real>()*(y*y*y - 4*z*z*z);
+ }
+ if (absx < 2)
+ {
+ Real y = 2 - absx;
+ return boost::math::constants::sixth<Real>()*y*y*y;
+ }
+ return (Real) 0;
+}
+
+template<class Real>
+Real b3_spline_prime(Real x)
+{
+ if (x < 0)
+ {
+ return -b3_spline_prime(-x);
+ }
+
+ if (x < 1)
+ {
+ return x*(3*boost::math::constants::half<Real>()*x - 2);
+ }
+ if (x < 2)
+ {
+ return -boost::math::constants::half<Real>()*(2 - x)*(2 - x);
+ }
+ return (Real) 0;
+}
+
+
+template <class Real>
+template <class BidiIterator>
+cubic_b_spline_imp<Real>::cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+ Real left_endpoint_derivative, Real right_endpoint_derivative) : m_a(left_endpoint), m_avg(0)
+{
+ using boost::math::constants::third;
+
+ std::size_t length = end_p - f;
+
+ if (length < 5)
+ {
+ if (boost::math::isnan(left_endpoint_derivative) || boost::math::isnan(right_endpoint_derivative))
+ {
+ throw std::logic_error("Interpolation using a cubic b spline with derivatives estimated at the endpoints requires at least 5 points.\n");
+ }
+ if (length < 3)
+ {
+ throw std::logic_error("Interpolation using a cubic b spline requires at least 3 points.\n");
+ }
+ }
+
+ if (boost::math::isnan(left_endpoint))
+ {
+ throw std::logic_error("Left endpoint is NAN; this is disallowed.\n");
+ }
+ if (left_endpoint + length*step_size >= (std::numeric_limits<Real>::max)())
+ {
+ throw std::logic_error("Right endpoint overflows the maximum representable number of the specified precision.\n");
+ }
+ if (step_size <= 0)
+ {
+ throw std::logic_error("The step size must be strictly > 0.\n");
+ }
+ // Storing the inverse of the stepsize does provide a measurable speedup.
+ // It's not huge, but nonetheless worthwhile.
+ m_h_inv = 1/step_size;
+
+ // Following Kress's notation, s'(a) = a1, s'(b) = b1
+ Real a1 = left_endpoint_derivative;
+ // See the finite-difference table on Wikipedia for reference on how
+ // to construct high-order estimates for one-sided derivatives:
+ // https://en.wikipedia.org/wiki/Finite_difference_coefficient#Forward_and_backward_finite_difference
+ // Here, we estimate then to O(h^4), as that is the maximum accuracy we could obtain from this method.
+ if (boost::math::isnan(a1))
+ {
+ // For simple functions (linear, quadratic, so on)
+ // almost all the error comes from derivative estimation.
+ // This does pairwise summation which gives us another digit of accuracy over naive summation.
+ Real t0 = 4*(f[1] + third<Real>()*f[3]);
+ Real t1 = -(25*third<Real>()*f[0] + f[4])/4 - 3*f[2];
+ a1 = m_h_inv*(t0 + t1);
+ }
+
+ Real b1 = right_endpoint_derivative;
+ if (boost::math::isnan(b1))
+ {
+ size_t n = length - 1;
+ Real t0 = 4*(f[n-3] + third<Real>()*f[n - 1]);
+ Real t1 = -(25*third<Real>()*f[n - 4] + f[n])/4 - 3*f[n - 2];
+
+ b1 = m_h_inv*(t0 + t1);
+ }
+
+ // s(x) = \sum \alpha_i B_{3}( (x- x_i - a)/h )
+ // Of course we must reindex from Kress's notation, since he uses negative indices which make C++ unhappy.
+ m_beta.resize(length + 2, std::numeric_limits<Real>::quiet_NaN());
+
+ // Since the splines have compact support, they decay to zero very fast outside the endpoints.
+ // This is often very annoying; we'd like to evaluate the interpolant a little bit outside the
+ // boundary [a,b] without massive error.
+ // A simple way to deal with this is just to subtract the DC component off the signal, so we need the average.
+ // This algorithm for computing the average is recommended in
+ // http://www.heikohoffmann.de/htmlthesis/node134.html
+ Real t = 1;
+ for (size_t i = 0; i < length; ++i)
+ {
+ if (boost::math::isnan(f[i]))
+ {
+ std::string err = "This function you are trying to interpolate is a nan at index " + std::to_string(i) + "\n";
+ throw std::logic_error(err);
+ }
+ m_avg += (f[i] - m_avg) / t;
+ t += 1;
+ }
+
+
+ // Now we must solve an almost-tridiagonal system, which requires O(N) operations.
+ // There are, in fact 5 diagonals, but they only differ from zero on the first and last row,
+ // so we can patch up the tridiagonal row reduction algorithm to deal with two special rows.
+ // See Kress, equations 8.41
+ // The the "tridiagonal" matrix is:
+ // 1 0 -1
+ // 1 4 1
+ // 1 4 1
+ // 1 4 1
+ // ....
+ // 1 4 1
+ // 1 0 -1
+ // Numerical estimate indicate that as N->Infinity, cond(A) -> 6.9, so this matrix is good.
+ std::vector<Real> rhs(length + 2, std::numeric_limits<Real>::quiet_NaN());
+ std::vector<Real> super_diagonal(length + 2, std::numeric_limits<Real>::quiet_NaN());
+
+ rhs[0] = -2*step_size*a1;
+ rhs[rhs.size() - 1] = -2*step_size*b1;
+
+ super_diagonal[0] = 0;
+
+ for(size_t i = 1; i < rhs.size() - 1; ++i)
+ {
+ rhs[i] = 6*(f[i - 1] - m_avg);
+ super_diagonal[i] = 1;
+ }
+
+
+ // One step of row reduction on the first row to patch up the 5-diagonal problem:
+ // 1 0 -1 | r0
+ // 1 4 1 | r1
+ // mapsto:
+ // 1 0 -1 | r0
+ // 0 4 2 | r1 - r0
+ // mapsto
+ // 1 0 -1 | r0
+ // 0 1 1/2| (r1 - r0)/4
+ super_diagonal[1] = 0.5;
+ rhs[1] = (rhs[1] - rhs[0])/4;
+
+ // Now do a tridiagonal row reduction the standard way, until just before the last row:
+ for (size_t i = 2; i < rhs.size() - 1; ++i)
+ {
+ Real diagonal = 4 - super_diagonal[i - 1];
+ rhs[i] = (rhs[i] - rhs[i - 1])/diagonal;
+ super_diagonal[i] /= diagonal;
+ }
+
+ // Now the last row, which is in the form
+ // 1 sd[n-3] 0 | rhs[n-3]
+ // 0 1 sd[n-2] | rhs[n-2]
+ // 1 0 -1 | rhs[n-1]
+ Real final_subdiag = -super_diagonal[rhs.size() - 3];
+ rhs[rhs.size() - 1] = (rhs[rhs.size() - 1] - rhs[rhs.size() - 3])/final_subdiag;
+ Real final_diag = -1/final_subdiag;
+ // Now we're here:
+ // 1 sd[n-3] 0 | rhs[n-3]
+ // 0 1 sd[n-2] | rhs[n-2]
+ // 0 1 final_diag | (rhs[n-1] - rhs[n-3])/diag
+
+ final_diag = final_diag - super_diagonal[rhs.size() - 2];
+ rhs[rhs.size() - 1] = rhs[rhs.size() - 1] - rhs[rhs.size() - 2];
+
+
+ // Back substitutions:
+ m_beta[rhs.size() - 1] = rhs[rhs.size() - 1]/final_diag;
+ for(size_t i = rhs.size() - 2; i > 0; --i)
+ {
+ m_beta[i] = rhs[i] - super_diagonal[i]*m_beta[i + 1];
+ }
+ m_beta[0] = m_beta[2] + rhs[0];
+}
+
+template<class Real>
+Real cubic_b_spline_imp<Real>::operator()(Real x) const
+{
+ // See Kress, 8.40: Since B3 has compact support, we don't have to sum over all terms,
+ // just the (at most 5) whose support overlaps the argument.
+ Real z = m_avg;
+ Real t = m_h_inv*(x - m_a) + 1;
+
+ using std::max;
+ using std::min;
+ using std::ceil;
+ using std::floor;
+
+ size_t k_min = (size_t) max(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2)));
+ size_t k_max = (size_t) max(min(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))), (long) 0);
+ for (size_t k = k_min; k <= k_max; ++k)
+ {
+ z += m_beta[k]*b3_spline(t - k);
+ }
+
+ return z;
+}
+
+template<class Real>
+Real cubic_b_spline_imp<Real>::prime(Real x) const
+{
+ Real z = 0;
+ Real t = m_h_inv*(x - m_a) + 1;
+
+ using std::max;
+ using std::min;
+ using std::ceil;
+ using std::floor;
+
+ size_t k_min = (size_t) max(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2)));
+ size_t k_max = (size_t) min(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2)));
+
+ for (size_t k = k_min; k <= k_max; ++k)
+ {
+ z += m_beta[k]*b3_spline_prime(t - k);
+ }
+ return z*m_h_inv;
+}
+
+}}}
+#endif