[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 {
namespace math {
namespace interpolators {
namespace detail {

template<class RandomAccessContainer>
class cubic_hermite_detail {
public:
    using Real = typename RandomAccessContainer::value_type;
    using Size = typename RandomAccessContainer::size_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;
    }

    Size 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;
    using Size = typename RandomAccessContainer::size_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_;
    }


    Size 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;
    using Size = typename RandomAccessContainer::size_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_;
    }


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