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 cubic_hermite_detail(RandomAccessContainer && x,RandomAccessContainer && y,RandomAccessContainer dydx)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 push_back(Real x,Real y,Real dydx)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 operator ()(Real x) const69 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+2*t) + s0*(x-x0)) 100 + t*t*(y1*(3-2*t) + dx*s1*(t-1)); 101 return y; 102 } 103 prime(Real x) const104 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 - s0*dx)/(dx*dx); 129 Real d2 = (s1 - s0)/(2*dx); 130 Real c2 = 3*d1 - 2*d2; 131 Real c3 = 2*(d2 - d1)/dx; 132 return s0 + 2*c2*(x-x0) + 3*c3*(x-x0)*(x-x0); 133 } 134 135 operator <<(std::ostream & os,const cubic_hermite_detail & m)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 size() const148 auto size() const 149 { 150 return x_.size(); 151 } 152 bytes() const153 int64_t bytes() const 154 { 155 return 3*x_.size()*sizeof(Real) + 3*sizeof(x_); 156 } 157 domain() const158 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 cardinal_cubic_hermite_detail(RandomAccessContainer && y,RandomAccessContainer dydx,Real x0,Real dx)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 operator ()(Real x) const201 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 unchecked_evaluation(Real x) const219 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 r*r*(y0*(1+2*t) + dy0*t) 233 + t*t*(y1*(3-2*t) - dy1*r); 234 } 235 prime(Real x) const236 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 unchecked_prime(Real x) const254 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 = 6*t*(1-t)*(y1 - y0) + (3*t*t - 4*t+1)*dy0 + t*(3*t-2)*dy1; 267 return dy*inv_dx_; 268 } 269 270 size() const271 auto size() const 272 { 273 return y_.size(); 274 } 275 bytes() const276 int64_t bytes() const 277 { 278 return 2*y_.size()*sizeof(Real) + 2*sizeof(y_) + 2*sizeof(Real); 279 } 280 domain() const281 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 cardinal_cubic_hermite_detail_aos(RandomAccessContainer && dat,Real x0,Real dx)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 operator ()(Real x) const324 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 unchecked_evaluation(Real x) const342 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 r*r*(y0*(1+2*t) + dy0*t) 363 + t*t*(y1*(3-2*t) - dy1*r); 364 } 365 prime(Real x) const366 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 unchecked_prime(Real x) const384 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 = 6*t*(1-t)*(y1 - y0) + (3*t*t - 4*t+1)*dy0 + t*(3*t-2)*dy1; 401 return dy*inv_dx_; 402 } 403 404 size() const405 auto size() const 406 { 407 return dat_.size(); 408 } 409 bytes() const410 int64_t bytes() const 411 { 412 return dat_.size()*dat_[0].size()*sizeof(Real) + sizeof(dat_) + 2*sizeof(Real); 413 } 414 domain() const415 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 431