[ceph.git] / ceph / src / boost / boost / math / interpolators / detail / cubic_hermite_detail.hpp

// Copyright Nick Thompson, 2020
// 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_DETAIL_CUBIC_HERMITE_DETAIL_HPP
#define BOOST_MATH_INTERPOLATORS_DETAIL_CUBIC_HERMITE_DETAIL_HPP
#include <stdexcept>
#include <algorithm>
#include <cmath>
#include <iostream>
#include <sstream>
#include <limits>

namespace boost::math::interpolators::detail {

template<class RandomAccessContainer>
class cubic_hermite_detail {
public:
    using Real = typename RandomAccessContainer::value_type;

    cubic_hermite_detail(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer dydx)
     : x_{std::move(x)}, y_{std::move(y)}, dydx_{std::move(dydx)}
    {
        using std::abs;
        using std::isnan;
        if (x_.size() != y_.size())
        {
            throw std::domain_error("There must be the same number of ordinates as abscissas.");
        }
        if (x_.size() != dydx_.size())
        {
            throw std::domain_error("There must be the same number of ordinates as derivative values.");
        }
        if (x_.size() < 2)
        {
            throw std::domain_error("Must be at least two data points.");
        }
        Real x0 = x_[0];
        for (size_t i = 1; i < x_.size(); ++i)
        {
            Real x1 = x_[i];
            if (x1 <= x0)
            {
                std::ostringstream oss;
                oss.precision(std::numeric_limits<Real>::digits10+3);
                oss << "Abscissas must be listed in strictly increasing order x0 < x1 < ... < x_{n-1}, ";
                oss << "but at x[" << i - 1 << "] = " << x0 << ", and x[" << i << "] = " << x1 << ".\n";
                throw std::domain_error(oss.str());
            }
            x0 = x1;
        }
    }

    void push_back(Real x, Real y, Real dydx)
    {
        using std::abs;
        using std::isnan;
        if (x <= x_.back())
        {
             throw std::domain_error("Calling push_back must preserve the monotonicity of the x's");
        }
        x_.push_back(x);
        y_.push_back(y);
        dydx_.push_back(dydx);
    }

    Real operator()(Real x) const
    {
        if  (x < x_[0] || x > x_.back())
        {
            std::ostringstream oss;
            oss.precision(std::numeric_limits<Real>::digits10+3);
            oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
                << x_[0] << ", " << x_.back() << "]";
            throw std::domain_error(oss.str());
        }
        // We need t := (x-x_k)/(x_{k+1}-x_k) \in [0,1) for this to work.
        // Sadly this neccessitates this loathesome check, otherwise we get t = 1 at x = xf.
        if (x == x_.back())
        {
            return y_.back();
        }

        auto it = std::upper_bound(x_.begin(), x_.end(), x);
        auto i = std::distance(x_.begin(), it) -1;
        Real x0 = *(it-1);
        Real x1 = *it;
        Real y0 = y_[i];
        Real y1 = y_[i+1];
        Real s0 = dydx_[i];
        Real s1 = dydx_[i+1];
        Real dx = (x1-x0);
        Real t = (x-x0)/dx;

        // See the section 'Representations' in the page
        // https://en.wikipedia.org/wiki/Cubic_Hermite_spline
        Real y = (1-t)*(1-t)*(y0*(1+2*t) + s0*(x-x0))
              + t*t*(y1*(3-2*t) + dx*s1*(t-1));
        return y;
    }

    Real prime(Real x) const
    {
        if  (x < x_[0] || x > x_.back())
        {
            std::ostringstream oss;
            oss.precision(std::numeric_limits<Real>::digits10+3);
            oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
                << x_[0] << ", " << x_.back() << "]";
            throw std::domain_error(oss.str());
        }
        if (x == x_.back())
        {
            return dydx_.back();
        }
        auto it = std::upper_bound(x_.begin(), x_.end(), x);
        auto i = std::distance(x_.begin(), it) -1;
        Real x0 = *(it-1);
        Real x1 = *it;
        Real y0 = y_[i];
        Real y1 = y_[i+1];
        Real s0 = dydx_[i];
        Real s1 = dydx_[i+1];
        Real dx = (x1-x0);

        Real d1 = (y1 - y0 - s0*dx)/(dx*dx);
        Real d2 = (s1 - s0)/(2*dx);
        Real c2 = 3*d1 - 2*d2;
        Real c3 = 2*(d2 - d1)/dx;
        return s0 + 2*c2*(x-x0) + 3*c3*(x-x0)*(x-x0); 
    }


    friend std::ostream& operator<<(std::ostream & os, const cubic_hermite_detail & m)
    {
        os << "(x,y,y') = {";
        for (size_t i = 0; i < m.x_.size() - 1; ++i)
        {
            os << "(" << m.x_[i] << ", " << m.y_[i] << ", " << m.dydx_[i] << "),  ";
        }
        auto n = m.x_.size()-1;
        os << "(" << m.x_[n] << ", " << m.y_[n] << ", " << m.dydx_[n] << ")}";
        return os;
    }

    auto size() const
    {
        return x_.size();
    }

    int64_t bytes() const
    {
        return 3*x_.size()*sizeof(Real) + 3*sizeof(x_);
    }

    std::pair<Real, Real> domain() const
    {
        return {x_.front(), x_.back()};
    }

    RandomAccessContainer x_;
    RandomAccessContainer y_;
    RandomAccessContainer dydx_;
};

template<class RandomAccessContainer>
class cardinal_cubic_hermite_detail {
public:
    using Real = typename RandomAccessContainer::value_type;

    cardinal_cubic_hermite_detail(RandomAccessContainer && y, RandomAccessContainer dydx, Real x0, Real dx)
    : y_{std::move(y)}, dy_{std::move(dydx)}, x0_{x0}, inv_dx_{1/dx}
    {
        using std::abs;
        using std::isnan;
        if (y_.size() != dy_.size())
        {
            throw std::domain_error("There must be the same number of derivatives as ordinates.");
        }
        if (y_.size() < 2)
        {
            throw std::domain_error("Must be at least two data points.");
        }
        if (dx <= 0)
        {
            throw std::domain_error("dx > 0 is required.");
        }

        for (auto & dy : dy_)
        {
            dy *= dx;
        }
    }

    // Why not implement push_back? It's awkward: If the buffer is circular, x0_ += dx_.
    // If the buffer is not circular, x0_ is unchanged.
    // We need a concept for circular_buffer!

    inline Real operator()(Real x) const
    {
        const Real xf = x0_ + (y_.size()-1)/inv_dx_;
        if  (x < x0_ || x > xf)
        {
            std::ostringstream oss;
            oss.precision(std::numeric_limits<Real>::digits10+3);
            oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
                << x0_ << ", " << xf << "]";
            throw std::domain_error(oss.str());
        }
        if (x == xf)
        {
            return y_.back();
        }
        return this->unchecked_evaluation(x);
    }

    inline Real unchecked_evaluation(Real x) const
    {
        using std::floor;
        Real s = (x-x0_)*inv_dx_;
        Real ii = floor(s);
        auto i = static_cast<decltype(y_.size())>(ii);
        Real t = s - ii;
        Real y0 = y_[i];
        Real y1 = y_[i+1];
        Real dy0 = dy_[i];
        Real dy1 = dy_[i+1];

        Real r = 1-t;
        return r*r*(y0*(1+2*t) + dy0*t)
              + t*t*(y1*(3-2*t) - dy1*r);
    }

    inline Real prime(Real x) const
    {
        const Real xf = x0_ + (y_.size()-1)/inv_dx_;
        if  (x < x0_ || x > xf)
        {
            std::ostringstream oss;
            oss.precision(std::numeric_limits<Real>::digits10+3);
            oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
                << x0_ << ", " << xf << "]";
            throw std::domain_error(oss.str());
        }
        if (x == xf)
        {
            return dy_.back()*inv_dx_;
        }
        return this->unchecked_prime(x);
    }

    inline Real unchecked_prime(Real x) const
    {
        using std::floor;
        Real s = (x-x0_)*inv_dx_;
        Real ii = floor(s);
        auto i = static_cast<decltype(y_.size())>(ii);
        Real t = s - ii;
        Real y0 = y_[i];
        Real y1 = y_[i+1];
        Real dy0 = dy_[i];
        Real dy1 = dy_[i+1];

        Real dy = 6*t*(1-t)*(y1 - y0)  + (3*t*t - 4*t+1)*dy0 + t*(3*t-2)*dy1;
        return dy*inv_dx_;
    }


    auto size() const
    {
        return y_.size();
    }

    int64_t bytes() const
    {
        return 2*y_.size()*sizeof(Real) + 2*sizeof(y_) + 2*sizeof(Real);
    }

    std::pair<Real, Real> domain() const
    {
        Real xf = x0_ + (y_.size()-1)/inv_dx_;
        return {x0_, xf};
    }

private:

    RandomAccessContainer y_;
    RandomAccessContainer dy_;
    Real x0_;
    Real inv_dx_;
};


template<class RandomAccessContainer>
class cardinal_cubic_hermite_detail_aos {
public:
    using Point = typename RandomAccessContainer::value_type;
    using Real = typename Point::value_type;

    cardinal_cubic_hermite_detail_aos(RandomAccessContainer && dat, Real x0, Real dx)
    : dat_{std::move(dat)}, x0_{x0}, inv_dx_{1/dx}
    {
        if (dat_.size() < 2)
        {
            throw std::domain_error("Must be at least two data points.");
        }
        if (dat_[0].size() != 2)
        {
            throw std::domain_error("Each datum must contain (y, y'), and nothing else.");
        }
        if (dx <= 0)
        {
            throw std::domain_error("dx > 0 is required.");
        }

        for (auto & d : dat_)
        {
            d[1] *= dx;
        }
    }

    inline Real operator()(Real x) const
    {
        const Real xf = x0_ + (dat_.size()-1)/inv_dx_;
        if  (x < x0_ || x > xf)
        {
            std::ostringstream oss;
            oss.precision(std::numeric_limits<Real>::digits10+3);
            oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
                << x0_ << ", " << xf << "]";
            throw std::domain_error(oss.str());
        }
        if (x == xf)
        {
            return dat_.back()[0];
        }
        return this->unchecked_evaluation(x);
    }

    inline Real unchecked_evaluation(Real x) const
    {
        using std::floor;
        Real s = (x-x0_)*inv_dx_;
        Real ii = floor(s);
        auto i = static_cast<decltype(dat_.size())>(ii);

        Real t = s - ii;
        // If we had infinite precision, this would never happen.
        // But we don't have infinite precision.
        if (t == 0)
        {
            return dat_[i][0];
        }
        Real y0 = dat_[i][0];
        Real y1 = dat_[i+1][0];
        Real dy0 = dat_[i][1];
        Real dy1 = dat_[i+1][1];

        Real r = 1-t;
        return r*r*(y0*(1+2*t) + dy0*t)
              + t*t*(y1*(3-2*t) - dy1*r);
    }

    inline Real prime(Real x) const
    {
        const Real xf = x0_ + (dat_.size()-1)/inv_dx_;
        if  (x < x0_ || x > xf)
        {
            std::ostringstream oss;
            oss.precision(std::numeric_limits<Real>::digits10+3);
            oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
                << x0_ << ", " << xf << "]";
            throw std::domain_error(oss.str());
        }
        if (x == xf)
        {
            return dat_.back()[1]*inv_dx_;
        }
        return this->unchecked_prime(x);
    }

    inline Real unchecked_prime(Real x) const
    {
        using std::floor;
        Real s = (x-x0_)*inv_dx_;
        Real ii = floor(s);
        auto i = static_cast<decltype(dat_.size())>(ii);
        Real t = s - ii;
        if (t == 0)
        {
            return dat_[i][1]*inv_dx_;
        }
        Real y0 = dat_[i][0];
        Real dy0 = dat_[i][1];
        Real y1 = dat_[i+1][0];
        Real dy1 = dat_[i+1][1];

        Real dy = 6*t*(1-t)*(y1 - y0)  + (3*t*t - 4*t+1)*dy0 + t*(3*t-2)*dy1;
        return dy*inv_dx_;
    }


    auto size() const
    {
        return dat_.size();
    }

    int64_t bytes() const
    {
        return dat_.size()*dat_[0].size()*sizeof(Real) + sizeof(dat_) + 2*sizeof(Real);
    }

    std::pair<Real, Real> domain() const
    {
        Real xf = x0_ + (dat_.size()-1)/inv_dx_;
        return {x0_, xf};
    }


private:
    RandomAccessContainer dat_;
    Real x0_;
    Real inv_dx_;
};


}
#endif
Commit	Line	Data
f67539c2 TL	1	// Copyright Nick Thompson, 2020
	2	// Use, modification and distribution are subject to the
	3	// Boost Software License, Version 1.0.
	4	// (See accompanying file LICENSE_1_0.txt
	5	// or copy at http://www.boost.org/LICENSE_1_0.txt)
	6
	7	#ifndef BOOST_MATH_INTERPOLATORS_DETAIL_CUBIC_HERMITE_DETAIL_HPP
	8	#define BOOST_MATH_INTERPOLATORS_DETAIL_CUBIC_HERMITE_DETAIL_HPP
	9	#include <stdexcept>
	10	#include <algorithm>
	11	#include <cmath>
	12	#include <iostream>
	13	#include <sstream>
	14	#include <limits>
	15
	16	namespace boost::math::interpolators::detail {
	17
	18	template<class RandomAccessContainer>
	19	class cubic_hermite_detail {
	20	public:
	21	using Real = typename RandomAccessContainer::value_type;
	22
	23	cubic_hermite_detail(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer dydx)
	24	: x_{std::move(x)}, y_{std::move(y)}, dydx_{std::move(dydx)}
	25	{
	26	using std::abs;
	27	using std::isnan;
	28	if (x_.size() != y_.size())
	29	{
	30	throw std::domain_error("There must be the same number of ordinates as abscissas.");
	31	}
	32	if (x_.size() != dydx_.size())
	33	{
	34	throw std::domain_error("There must be the same number of ordinates as derivative values.");
	35	}
	36	if (x_.size() < 2)
	37	{
	38	throw std::domain_error("Must be at least two data points.");
	39	}
	40	Real x0 = x_[0];
	41	for (size_t i = 1; i < x_.size(); ++i)
	42	{
	43	Real x1 = x_[i];
	44	if (x1 <= x0)
	45	{
	46	std::ostringstream oss;
	47	oss.precision(std::numeric_limits<Real>::digits10+3);
	48	oss << "Abscissas must be listed in strictly increasing order x0 < x1 < ... < x_{n-1}, ";
	49	oss << "but at x[" << i - 1 << "] = " << x0 << ", and x[" << i << "] = " << x1 << ".\n";
	50	throw std::domain_error(oss.str());
	51	}
	52	x0 = x1;
	53	}
	54	}
	55
	56	void push_back(Real x, Real y, Real dydx)
	57	{
	58	using std::abs;
	59	using std::isnan;
	60	if (x <= x_.back())
	61	{
	62	throw std::domain_error("Calling push_back must preserve the monotonicity of the x's");
	63	}
	64	x_.push_back(x);
65	y_.push_back(y);
66	dydx_.push_back(dydx);
67	}
68
69	Real operator()(Real x) const
70	{
71	if (x < x_[0] \|\| x > x_.back())
72	{
73	std::ostringstream oss;
74	oss.precision(std::numeric_limits<Real>::digits10+3);
75	oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
76	<< x_[0] << ", " << x_.back() << "]";
77	throw std::domain_error(oss.str());
78	}
79	// We need t := (x-x_k)/(x_{k+1}-x_k) \in [0,1) for this to work.
80	// Sadly this neccessitates this loathesome check, otherwise we get t = 1 at x = xf.
81	if (x == x_.back())
82	{
83	return y_.back();
84	}
85
86	auto it = std::upper_bound(x_.begin(), x_.end(), x);
87	auto i = std::distance(x_.begin(), it) -1;
88	Real x0 = *(it-1);
89	Real x1 = *it;
90	Real y0 = y_[i];
91	Real y1 = y_[i+1];
92	Real s0 = dydx_[i];
93	Real s1 = dydx_[i+1];
94	Real dx = (x1-x0);
95	Real t = (x-x0)/dx;
96
97	// See the section 'Representations' in the page
98	// https://en.wikipedia.org/wiki/Cubic_Hermite_spline
99	Real y = (1-t)(1-t)(y0(1+2t) + s0*(x-x0))
100	+ tt(y1(3-2t) + dxs1(t-1));
101	return y;
102	}
103
104	Real prime(Real x) const
105	{
106	if (x < x_[0] \|\| x > x_.back())
107	{
108	std::ostringstream oss;
109	oss.precision(std::numeric_limits<Real>::digits10+3);
110	oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
111	<< x_[0] << ", " << x_.back() << "]";
112	throw std::domain_error(oss.str());
113	}
114	if (x == x_.back())
115	{
116	return dydx_.back();
117	}
118	auto it = std::upper_bound(x_.begin(), x_.end(), x);
119	auto i = std::distance(x_.begin(), it) -1;
120	Real x0 = *(it-1);
121	Real x1 = *it;
122	Real y0 = y_[i];
123	Real y1 = y_[i+1];
124	Real s0 = dydx_[i];
125	Real s1 = dydx_[i+1];
126	Real dx = (x1-x0);
127
128	Real d1 = (y1 - y0 - s0dx)/(dxdx);
129	Real d2 = (s1 - s0)/(2*dx);
130	Real c2 = 3d1 - 2d2;
131	Real c3 = 2*(d2 - d1)/dx;
132	return s0 + 2c2(x-x0) + 3c3(x-x0)*(x-x0);
133	}
134
135
136	friend std::ostream& operator<<(std::ostream & os, const cubic_hermite_detail & m)
137	{
138	os << "(x,y,y') = {";
139	for (size_t i = 0; i < m.x_.size() - 1; ++i)
140	{
141	os << "(" << m.x_[i] << ", " << m.y_[i] << ", " << m.dydx_[i] << "), ";
142	}
143	auto n = m.x_.size()-1;
144	os << "(" << m.x_[n] << ", " << m.y_[n] << ", " << m.dydx_[n] << ")}";
145	return os;
146	}
147
148	auto size() const
149	{
150	return x_.size();
151	}
152
153	int64_t bytes() const
154	{
155	return 3x_.size()sizeof(Real) + 3*sizeof(x_);
156	}
157
158	std::pair<Real, Real> domain() const
159	{
160	return {x_.front(), x_.back()};
161	}
162
163	RandomAccessContainer x_;
164	RandomAccessContainer y_;
165	RandomAccessContainer dydx_;
166	};
167
168	template<class RandomAccessContainer>
169	class cardinal_cubic_hermite_detail {
170	public:
171	using Real = typename RandomAccessContainer::value_type;
172
173	cardinal_cubic_hermite_detail(RandomAccessContainer && y, RandomAccessContainer dydx, Real x0, Real dx)
174	: y_{std::move(y)}, dy_{std::move(dydx)}, x0_{x0}, inv_dx_{1/dx}
175	{
176	using std::abs;
177	using std::isnan;
178	if (y_.size() != dy_.size())
179	{
180	throw std::domain_error("There must be the same number of derivatives as ordinates.");
181	}
182	if (y_.size() < 2)
183	{
184	throw std::domain_error("Must be at least two data points.");
185	}
186	if (dx <= 0)
187	{
188	throw std::domain_error("dx > 0 is required.");
189	}
190
191	for (auto & dy : dy_)
192	{
193	dy *= dx;
194	}
195	}
196
197	// Why not implement push_back? It's awkward: If the buffer is circular, x0_ += dx_.
198	// If the buffer is not circular, x0_ is unchanged.
199	// We need a concept for circular_buffer!
200
201	inline Real operator()(Real x) const
202	{
203	const Real xf = x0_ + (y_.size()-1)/inv_dx_;
204	if (x < x0_ \|\| x > xf)
205	{
206	std::ostringstream oss;
207	oss.precision(std::numeric_limits<Real>::digits10+3);
208	oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
209	<< x0_ << ", " << xf << "]";
210	throw std::domain_error(oss.str());
211	}
212	if (x == xf)
213	{
214	return y_.back();
215	}
216	return this->unchecked_evaluation(x);
217	}
218
219	inline Real unchecked_evaluation(Real x) const
220	{
221	using std::floor;
222	Real s = (x-x0_)*inv_dx_;
223	Real ii = floor(s);
224	auto i = static_cast<decltype(y_.size())>(ii);
225	Real t = s - ii;
226	Real y0 = y_[i];
227	Real y1 = y_[i+1];
228	Real dy0 = dy_[i];
229	Real dy1 = dy_[i+1];
230
231	Real r = 1-t;
232	return rr(y0(1+2t) + dy0*t)
233	+ tt(y1(3-2t) - dy1*r);
234	}
235
236	inline Real prime(Real x) const
237	{
238	const Real xf = x0_ + (y_.size()-1)/inv_dx_;
239	if (x < x0_ \|\| x > xf)
240	{
241	std::ostringstream oss;
242	oss.precision(std::numeric_limits<Real>::digits10+3);
243	oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
244	<< x0_ << ", " << xf << "]";
245	throw std::domain_error(oss.str());
246	}
247	if (x == xf)
248	{
249	return dy_.back()*inv_dx_;
250	}
251	return this->unchecked_prime(x);
252	}
253
254	inline Real unchecked_prime(Real x) const
255	{
256	using std::floor;
257	Real s = (x-x0_)*inv_dx_;
258	Real ii = floor(s);
259	auto i = static_cast<decltype(y_.size())>(ii);
260	Real t = s - ii;
261	Real y0 = y_[i];
262	Real y1 = y_[i+1];
263	Real dy0 = dy_[i];
264	Real dy1 = dy_[i+1];
265
266	Real dy = 6t(1-t)(y1 - y0) + (3tt - 4t+1)dy0 + t(3t-2)dy1;
267	return dy*inv_dx_;
268	}
269
270
271	auto size() const
272	{
273	return y_.size();
274	}
275
276	int64_t bytes() const
277	{
278	return 2y_.size()sizeof(Real) + 2sizeof(y_) + 2sizeof(Real);
279	}
280
281	std::pair<Real, Real> domain() const
282	{
283	Real xf = x0_ + (y_.size()-1)/inv_dx_;
284	return {x0_, xf};
285	}
286
287	private:
288
289	RandomAccessContainer y_;
290	RandomAccessContainer dy_;
291	Real x0_;
292	Real inv_dx_;
293	};
294
295
296	template<class RandomAccessContainer>
297	class cardinal_cubic_hermite_detail_aos {
298	public:
299	using Point = typename RandomAccessContainer::value_type;
300	using Real = typename Point::value_type;
301
302	cardinal_cubic_hermite_detail_aos(RandomAccessContainer && dat, Real x0, Real dx)
303	: dat_{std::move(dat)}, x0_{x0}, inv_dx_{1/dx}
304	{
305	if (dat_.size() < 2)
306	{
307	throw std::domain_error("Must be at least two data points.");
308	}
309	if (dat_[0].size() != 2)
310	{
311	throw std::domain_error("Each datum must contain (y, y'), and nothing else.");
312	}
313	if (dx <= 0)
314	{
315	throw std::domain_error("dx > 0 is required.");
316	}
317
318	for (auto & d : dat_)
319	{
320	d[1] *= dx;
321	}
322	}
323
324	inline Real operator()(Real x) const
325	{
326	const Real xf = x0_ + (dat_.size()-1)/inv_dx_;
327	if (x < x0_ \|\| x > xf)
328	{
329	std::ostringstream oss;
330	oss.precision(std::numeric_limits<Real>::digits10+3);
331	oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
332	<< x0_ << ", " << xf << "]";
333	throw std::domain_error(oss.str());
334	}
335	if (x == xf)
336	{
337	return dat_.back()[0];
338	}
339	return this->unchecked_evaluation(x);
340	}
341
342	inline Real unchecked_evaluation(Real x) const
343	{
344	using std::floor;
345	Real s = (x-x0_)*inv_dx_;
346	Real ii = floor(s);
347	auto i = static_cast<decltype(dat_.size())>(ii);
348
349	Real t = s - ii;
350	// If we had infinite precision, this would never happen.
351	// But we don't have infinite precision.
352	if (t == 0)
353	{
354	return dat_[i][0];
355	}
356	Real y0 = dat_[i][0];
357	Real y1 = dat_[i+1][0];
358	Real dy0 = dat_[i][1];
359	Real dy1 = dat_[i+1][1];
360
361	Real r = 1-t;
362	return rr(y0(1+2t) + dy0*t)
363	+ tt(y1(3-2t) - dy1*r);
364	}
365
366	inline Real prime(Real x) const
367	{
368	const Real xf = x0_ + (dat_.size()-1)/inv_dx_;
369	if (x < x0_ \|\| x > xf)
370	{
371	std::ostringstream oss;
372	oss.precision(std::numeric_limits<Real>::digits10+3);
373	oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
374	<< x0_ << ", " << xf << "]";
375	throw std::domain_error(oss.str());
376	}
377	if (x == xf)
378	{
379	return dat_.back()[1]*inv_dx_;
380	}
381	return this->unchecked_prime(x);
382	}
383
384	inline Real unchecked_prime(Real x) const
385	{
386	using std::floor;
387	Real s = (x-x0_)*inv_dx_;
388	Real ii = floor(s);
389	auto i = static_cast<decltype(dat_.size())>(ii);
390	Real t = s - ii;
391	if (t == 0)
392	{
393	return dat_[i][1]*inv_dx_;
394	}
395	Real y0 = dat_[i][0];
396	Real dy0 = dat_[i][1];
397	Real y1 = dat_[i+1][0];
398	Real dy1 = dat_[i+1][1];
399
400	Real dy = 6t(1-t)(y1 - y0) + (3tt - 4t+1)dy0 + t(3t-2)dy1;
401	return dy*inv_dx_;
402	}
403
404
405	auto size() const
406	{
407	return dat_.size();
408	}
409
410	int64_t bytes() const
411	{
412	return dat_.size()dat_[0].size()sizeof(Real) + sizeof(dat_) + 2*sizeof(Real);
413	}
414
415	std::pair<Real, Real> domain() const
416	{
417	Real xf = x0_ + (dat_.size()-1)/inv_dx_;
418	return {x0_, xf};
419	}
420
421
422	private:
423	RandomAccessContainer dat_;
424	Real x0_;
425	Real inv_dx_;
426	};
427
428
429	}
430	#endif