interpolators/detail/whittaker_shannon_detail.hpp

0001 // Copyright Nick Thompson, 2019
0002 // Use, modification and distribution are subject to the
0003 // Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006 #ifndef BOOST_MATH_INTERPOLATORS_WHITAKKER_SHANNON_DETAIL_HPP
0007 #define BOOST_MATH_INTERPOLATORS_WHITAKKER_SHANNON_DETAIL_HPP
0008 #include <cmath>
0009 #include <boost/math/tools/assert.hpp>
0010 #include <boost/math/constants/constants.hpp>
0011 #include <boost/math/special_functions/sin_pi.hpp>
0012 #include <boost/math/special_functions/cos_pi.hpp>
0013
0014 namespace boost { namespace math { namespace interpolators { namespace detail {
0015
0016 template<class RandomAccessContainer>
0017 class whittaker_shannon_detail {
0018 public:
0019
0020     using Real = typename RandomAccessContainer::value_type;
0021     whittaker_shannon_detail(RandomAccessContainer&& y, Real const & t0, Real const & h) : m_y{std::move(y)}, m_t0{t0}, m_h{h}
0022     {
0023         for (size_t i = 1; i < m_y.size(); i += 2)
0024         {
0025             m_y[i] = -m_y[i];
0026         }
0027     }
0028
0029     inline Real operator()(Real t) const {
0030         using boost::math::constants::pi;
0031         using std::isfinite;
0032         using std::floor;
0033         using std::ceil;
0034         Real y = 0;
0035         Real x = (t - m_t0)/m_h;
0036         Real z = x;
0037         auto it = m_y.begin();
0038
0039         // For some reason, neither clang nor g++ will cache the address of m_y.end() in a register.
0040         // Hence make a copy of it:
0041         auto end = m_y.end();
0042         while(it != end)
0043         {
0044             y += *it++/z;
0045             z -= 1;
0046         }
0047
0048         if (!isfinite(y))
0049         {
0050             BOOST_MATH_ASSERT_MSG(floor(x) == ceil(x), "Floor and ceiling should be equal.\n");
0051             auto i = static_cast<size_t>(floor(x));
0052             if (i & 1)
0053             {
0054                 return -m_y[i];
0055             }
0056             return m_y[i];
0057         }
0058         return y*boost::math::sin_pi(x)/pi<Real>();
0059     }
0060
0061     Real prime(Real t) const {
0062         using boost::math::constants::pi;
0063         using std::isfinite;
0064         using std::floor;
0065         using std::ceil;
0066
0067         Real x = (t - m_t0)/m_h;
0068         if (ceil(x) == x) {
0069             Real s = 0;
0070             auto j = static_cast<long>(x);
0071             auto n = static_cast<long>(m_y.size());
0072             for (long i = 0; i < n; ++i)
0073             {
0074                 if (j - i != 0)
0075                 {
0076                     s += m_y[i]/(j-i);
0077                 }
0078                 // else derivative of sinc at zero is zero.
0079             }
0080             if (j & 1) {
0081                 s /= -m_h;
0082             } else {
0083                 s /= m_h;
0084             }
0085             return s;
0086         }
0087         Real z = x;
0088         auto it = m_y.begin();
0089         Real cospix = boost::math::cos_pi(x);
0090         Real sinpix_div_pi = boost::math::sin_pi(x)/pi<Real>();
0091
0092         Real s = 0;
0093         auto end = m_y.end();
0094         while(it != end)
0095         {
0096             s += (*it++)*(z*cospix - sinpix_div_pi)/(z*z);
0097             z -= 1;
0098         }
0099
0100         return s/m_h;
0101     }
0102
0103
0104
0105     Real operator[](size_t i) const {
0106         if (i & 1)
0107         {
0108             return -m_y[i];
0109         }
0110         return m_y[i];
0111     }
0112
0113     RandomAccessContainer&& return_data() {
0114         for (size_t i = 1; i < m_y.size(); i += 2)
0115         {
0116             m_y[i] = -m_y[i];
0117         }
0118         return std::move(m_y);
0119     }
0120
0121
0122 private:
0123     RandomAccessContainer m_y;
0124     Real m_t0;
0125     Real m_h;
0126 };
0127 }}}}
0128 #endif