ceph/src/boost/boost/math/interpolators/detail/quintic_hermite_detail.hpp

   1 /*
   2  * Copyright Nick Thompson, 2020
   3  * Use, modification and distribution are subject to the
   4  * Boost Software License, Version 1.0. (See accompanying file
   5  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6  */
   7 #ifndef BOOST_MATH_INTERPOLATORS_DETAIL_QUINTIC_HERMITE_DETAIL_HPP
   8 #define BOOST_MATH_INTERPOLATORS_DETAIL_QUINTIC_HERMITE_DETAIL_HPP
   9 #include <algorithm>
  10 #include <stdexcept>
  11 #include <sstream>
  12 #include <cmath>
  13
  14 namespace boost {
  15 namespace math {
  16 namespace interpolators {
  17 namespace detail {
  18
  19 template<class RandomAccessContainer>
  20 class quintic_hermite_detail {
  21 public:
  22     using Real = typename RandomAccessContainer::value_type;
  23     quintic_hermite_detail(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2) : x_{std::move(x)}, y_{std::move(y)}, dydx_{std::move(dydx)}, d2ydx2_{std::move(d2ydx2)}
  24     {
  25         if (x_.size() != y_.size())
  26         {
  27             throw std::domain_error("Number of abscissas must = number of ordinates.");
  28         }
  29         if (x_.size() != dydx_.size())
  30         {
  31             throw std::domain_error("Numbers of derivatives must = number of abscissas.");
  32         }
  33         if (x_.size() != d2ydx2_.size())
  34         {
  35             throw std::domain_error("Number of second derivatives must equal number of abscissas.");
  36         }
  37         if (x_.size() < 2)
  38         {
  39             throw std::domain_error("At least 2 abscissas are required.");
  40         }
  41         Real x0 = x_[0];
  42         for (decltype(x_.size()) i = 1; i < x_.size(); ++i)
  43         {
  44             Real x1 = x_[i];
  45             if (x1 <= x0)
  46             {
  47                 throw std::domain_error("Abscissas must be sorted in strictly increasing order x0 < x1 < ... < x_{n-1}");
  48             }
  49             x0 = x1;
  50         }
  51     }
  52
  53     void push_back(Real x, Real y, Real dydx, Real d2ydx2)
  54     {
  55         using std::abs;
  56         using std::isnan;
  57         if (x <= x_.back())
  58         {
  59              throw std::domain_error("Calling push_back must preserve the monotonicity of the x's");
  60         }
  61         x_.push_back(x);
  62         y_.push_back(y);
  63         dydx_.push_back(dydx);
  64         d2ydx2_.push_back(d2ydx2);
  65     }
  66
  67     inline Real operator()(Real x) const
  68     {
  69         if  (x < x_[0] || x > x_.back())
  70         {
  71             std::ostringstream oss;
  72             oss.precision(std::numeric_limits<Real>::digits10+3);
  73             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
  74                 << x_[0] << ", " << x_.back() << "]";
  75             throw std::domain_error(oss.str());
  76         }
  77         // We need t := (x-x_k)/(x_{k+1}-x_k) \in [0,1) for this to work.
  78         // Sadly this neccessitates this loathesome check, otherwise we get t = 1 at x = xf.
  79         if (x == x_.back())
  80         {
  81             return y_.back();
  82         }
  83
  84         auto it = std::upper_bound(x_.begin(), x_.end(), x);
  85         auto i = std::distance(x_.begin(), it) -1;
  86         Real x0 = *(it-1);
  87         Real x1 = *it;
  88         Real y0 = y_[i];
  89         Real y1 = y_[i+1];
  90         Real v0 = dydx_[i];
  91         Real v1 = dydx_[i+1];
  92         Real a0 = d2ydx2_[i];
  93         Real a1 = d2ydx2_[i+1];
  94
  95         Real dx = (x1-x0);
  96         Real t = (x-x0)/dx;
  97         Real t2 = t*t;
  98         Real t3 = t2*t;
  99
 100         // See the 'Basis functions' section of:
 101         // https://www.rose-hulman.edu/~finn/CCLI/Notes/day09.pdf
 102         // Also: https://github.com/MrHexxx/QuinticHermiteSpline/blob/master/HermiteSpline.cs
 103         Real y = (1- t3*(10 + t*(-15 + 6*t)))*y0;
 104         y += t*(1+ t2*(-6 + t*(8 -3*t)))*v0*dx;
 105         y += t2*(1 + t*(-3 + t*(3-t)))*a0*dx*dx/2;
 106         y += t3*((1 + t*(-2 + t))*a1*dx*dx/2 + (-4 + t*(7 - 3*t))*v1*dx + (10 + t*(-15 + 6*t))*y1);
 107         return y;
 108     }
 109
 110     inline Real prime(Real x) const
 111     {
 112         if  (x < x_[0] || x > x_.back())
 113         {
 114             std::ostringstream oss;
 115             oss.precision(std::numeric_limits<Real>::digits10+3);
 116             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
 117                 << x_[0] << ", " << x_.back() << "]";
 118             throw std::domain_error(oss.str());
 119         }
 120         if (x == x_.back())
 121         {
 122             return dydx_.back();
 123         }
 124
 125         auto it = std::upper_bound(x_.begin(), x_.end(), x);
 126         auto i = std::distance(x_.begin(), it) -1;
 127         Real x0 = *(it-1);
 128         Real x1 = *it;
 129         Real dx = x1 - x0;
 130
 131         Real y0 = y_[i];
 132         Real y1 = y_[i+1];
 133         Real v0 = dydx_[i];
 134         Real v1 = dydx_[i+1];
 135         Real a0 = d2ydx2_[i];
 136         Real a1 = d2ydx2_[i+1];
 137         Real t= (x-x0)/dx;
 138         Real t2 = t*t;
 139
 140         Real dydx = 30*t2*(1 - 2*t + t*t)*(y1-y0)/dx;
 141         dydx += (1-18*t*t + 32*t*t*t - 15*t*t*t*t)*v0 - t*t*(12 - 28*t + 15*t*t)*v1;
 142         dydx += (t*dx/2)*((2 - 9*t + 12*t*t - 5*t*t*t)*a0 + t*(3 - 8*t + 5*t*t)*a1);
 143         return dydx;
 144     }
 145
 146     inline Real double_prime(Real x) const
 147     {
 148         if  (x < x_[0] || x > x_.back())
 149         {
 150             std::ostringstream oss;
 151             oss.precision(std::numeric_limits<Real>::digits10+3);
 152             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
 153                 << x_[0] << ", " << x_.back() << "]";
 154             throw std::domain_error(oss.str());
 155         }
 156         if (x == x_.back())
 157         {
 158             return d2ydx2_.back();
 159         }
 160
 161         auto it = std::upper_bound(x_.begin(), x_.end(), x);
 162         auto i = std::distance(x_.begin(), it) -1;
 163         Real x0 = *(it-1);
 164         Real x1 = *it;
 165         Real dx = x1 - x0;
 166
 167         Real y0 = y_[i];
 168         Real y1 = y_[i+1];
 169         Real v0 = dydx_[i];
 170         Real v1 = dydx_[i+1];
 171         Real a0 = d2ydx2_[i];
 172         Real a1 = d2ydx2_[i+1];
 173         Real t = (x-x0)/dx;
 174
 175         Real d2ydx2 = 60*t*(1 + t*(-3 + 2*t))*(y1-y0)/(dx*dx);
 176         d2ydx2 += 12*t*(-3 + t*(8 - 5*t))*v0/dx;
 177         d2ydx2 -= 12*t*(2 + t*(-7 + 5*t))*v1/dx;
 178         d2ydx2 += (1 + t*(-9 + t*(18 - 10*t)))*a0;
 179         d2ydx2 += t*(3 + t*(-12 + 10*t))*a1;
 180
 181         return d2ydx2;
 182     }
 183
 184     friend std::ostream& operator<<(std::ostream & os, const quintic_hermite_detail & m)
 185     {
 186         os << "(x,y,y') = {";
 187         for (size_t i = 0; i < m.x_.size() - 1; ++i) {
 188             os << "(" << m.x_[i] << ", " << m.y_[i] << ", " << m.dydx_[i] << ", " << m.d2ydx2_[i] << "),  ";
 189         }
 190         auto n = m.x_.size()-1;
 191         os << "(" << m.x_[n] << ", " << m.y_[n] << ", " << m.dydx_[n] << ", " << m.d2ydx2_[n] << ")}";
 192         return os;
 193     }
 194
 195     int64_t bytes() const
 196     {
 197         return 4*x_.size()*sizeof(x_);
 198     }
 199
 200     std::pair<Real, Real> domain() const
 201     {
 202         return {x_.front(), x_.back()};
 203     }
 204
 205 private:
 206     RandomAccessContainer x_;
 207     RandomAccessContainer y_;
 208     RandomAccessContainer dydx_;
 209     RandomAccessContainer d2ydx2_;
 210 };
 211
 212
 213 template<class RandomAccessContainer>
 214 class cardinal_quintic_hermite_detail {
 215 public:
 216     using Real = typename RandomAccessContainer::value_type;
 217     cardinal_quintic_hermite_detail(RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2, Real x0, Real dx)
 218     : y_{std::move(y)}, dy_{std::move(dydx)}, d2y_{std::move(d2ydx2)}, x0_{x0}, inv_dx_{1/dx}
 219     {
 220         if (y_.size() != dy_.size())
 221         {
 222             throw std::domain_error("Numbers of derivatives must = number of abscissas.");
 223         }
 224         if (y_.size() != d2y_.size())
 225         {
 226             throw std::domain_error("Number of second derivatives must equal number of abscissas.");
 227         }
 228         if (y_.size() < 2)
 229         {
 230             throw std::domain_error("At least 2 abscissas are required.");
 231         }
 232         if (dx <= 0)
 233         {
 234             throw std::domain_error("dx > 0 is required.");
 235         }
 236
 237         for (auto & dy : dy_)
 238         {
 239             dy *= dx;
 240         }
 241
 242         for (auto & d2y : d2y_)
 243         {
 244             d2y *= (dx*dx)/2;
 245         }
 246     }
 247
 248
 249     inline Real operator()(Real x) const
 250     {
 251         const Real xf = x0_ + (y_.size()-1)/inv_dx_;
 252         if  (x < x0_ || x > xf)
 253         {
 254             std::ostringstream oss;
 255             oss.precision(std::numeric_limits<Real>::digits10+3);
 256             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
 257                 << x0_ << ", " << xf << "]";
 258             throw std::domain_error(oss.str());
 259         }
 260         if (x == xf)
 261         {
 262             return y_.back();
 263         }
 264         return this->unchecked_evaluation(x);
 265     }
 266
 267     inline Real unchecked_evaluation(Real x) const
 268     {
 269         using std::floor;
 270         Real s = (x-x0_)*inv_dx_;
 271         Real ii = floor(s);
 272         auto i = static_cast<decltype(y_.size())>(ii);
 273         Real t = s - ii;
 274         if (t == 0)
 275         {
 276             return y_[i];
 277         }
 278         Real y0 = y_[i];
 279         Real y1 = y_[i+1];
 280         Real dy0 = dy_[i];
 281         Real dy1 = dy_[i+1];
 282         Real d2y0 = d2y_[i];
 283         Real d2y1 = d2y_[i+1];
 284
 285         // See the 'Basis functions' section of:
 286         // https://www.rose-hulman.edu/~finn/CCLI/Notes/day09.pdf
 287         // Also: https://github.com/MrHexxx/QuinticHermiteSpline/blob/master/HermiteSpline.cs
 288         Real y = (1- t*t*t*(10 + t*(-15 + 6*t)))*y0;
 289         y += t*(1+ t*t*(-6 + t*(8 -3*t)))*dy0;
 290         y += t*t*(1 + t*(-3 + t*(3-t)))*d2y0;
 291         y += t*t*t*((1 + t*(-2 + t))*d2y1 + (-4 + t*(7 -3*t))*dy1 + (10 + t*(-15 + 6*t))*y1);
 292         return y;
 293     }
 294
 295     inline Real prime(Real x) const
 296     {
 297         const Real xf = x0_ + (y_.size()-1)/inv_dx_;
 298         if  (x < x0_ || x > xf)
 299         {
 300             std::ostringstream oss;
 301             oss.precision(std::numeric_limits<Real>::digits10+3);
 302             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
 303                 << x0_ << ", " << xf << "]";
 304             throw std::domain_error(oss.str());
 305         }
 306         if (x == xf)
 307         {
 308             return dy_.back()*inv_dx_;
 309         }
 310
 311         return this->unchecked_prime(x);
 312     }
 313
 314     inline Real unchecked_prime(Real x) const
 315     {
 316         using std::floor;
 317         Real s = (x-x0_)*inv_dx_;
 318         Real ii = floor(s);
 319         auto i = static_cast<decltype(y_.size())>(ii);
 320         Real t = s - ii;
 321         if (t == 0)
 322         {
 323             return dy_[i]*inv_dx_;
 324         }
 325         Real y0 = y_[i];
 326         Real y1 = y_[i+1];
 327         Real dy0 = dy_[i];
 328         Real dy1 = dy_[i+1];
 329         Real d2y0 = d2y_[i];
 330         Real d2y1 = d2y_[i+1];
 331
 332         Real dydx = 30*t*t*(1 - 2*t + t*t)*(y1-y0);
 333         dydx += (1-18*t*t + 32*t*t*t - 15*t*t*t*t)*dy0 - t*t*(12 - 28*t + 15*t*t)*dy1;
 334         dydx += t*((2 - 9*t + 12*t*t - 5*t*t*t)*d2y0 + t*(3 - 8*t + 5*t*t)*d2y1);
 335         dydx *= inv_dx_;
 336         return dydx;
 337     }
 338
 339     inline Real double_prime(Real x) const
 340     {
 341         const Real xf = x0_ + (y_.size()-1)/inv_dx_;
 342         if  (x < x0_ || x > xf) {
 343             std::ostringstream oss;
 344             oss.precision(std::numeric_limits<Real>::digits10+3);
 345             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
 346                 << x0_ << ", " << xf << "]";
 347             throw std::domain_error(oss.str());
 348         }
 349         if (x == xf)
 350         {
 351             return d2y_.back()*2*inv_dx_*inv_dx_;
 352         }
 353
 354         return this->unchecked_double_prime(x);
 355     }
 356
 357     inline Real unchecked_double_prime(Real x) const
 358     {
 359         using std::floor;
 360         Real s = (x-x0_)*inv_dx_;
 361         Real ii = floor(s);
 362         auto i = static_cast<decltype(y_.size())>(ii);
 363         Real t = s - ii;
 364         if (t==0)
 365         {
 366             return d2y_[i]*2*inv_dx_*inv_dx_;
 367         }
 368
 369         Real y0 = y_[i];
 370         Real y1 = y_[i+1];
 371         Real dy0 = dy_[i];
 372         Real dy1 = dy_[i+1];
 373         Real d2y0 = d2y_[i];
 374         Real d2y1 = d2y_[i+1];
 375
 376         Real d2ydx2 = 60*t*(1 - 3*t + 2*t*t)*(y1 - y0)*inv_dx_*inv_dx_;
 377         d2ydx2 += (12*t)*((-3 + 8*t - 5*t*t)*dy0 - (2 - 7*t + 5*t*t)*dy1);
 378         d2ydx2 += (1 - 9*t + 18*t*t - 10*t*t*t)*d2y0*(2*inv_dx_*inv_dx_) + t*(3 - 12*t + 10*t*t)*d2y1*(2*inv_dx_*inv_dx_);
 379         return d2ydx2;
 380     }
 381
 382     int64_t bytes() const
 383     {
 384         return 3*y_.size()*sizeof(Real) + 2*sizeof(Real);
 385     }
 386
 387     std::pair<Real, Real> domain() const
 388     {
 389         Real xf = x0_ + (y_.size()-1)/inv_dx_;
 390         return {x0_, xf};
 391     }
 392
 393 private:
 394     RandomAccessContainer y_;
 395     RandomAccessContainer dy_;
 396     RandomAccessContainer d2y_;
 397     Real x0_;
 398     Real inv_dx_;
 399 };
 400
 401
 402 template<class RandomAccessContainer>
 403 class cardinal_quintic_hermite_detail_aos {
 404 public:
 405     using Point = typename RandomAccessContainer::value_type;
 406     using Real = typename Point::value_type;
 407     cardinal_quintic_hermite_detail_aos(RandomAccessContainer && data, Real x0, Real dx)
 408     : data_{std::move(data)} , x0_{x0}, inv_dx_{1/dx}
 409     {
 410         if (data_.size() < 2)
 411         {
 412             throw std::domain_error("At least two points are required to interpolate using cardinal_quintic_hermite_aos");
 413         }
 414
 415         if (data_[0].size() != 3)
 416         {
 417             throw std::domain_error("Each datum passed to the cardinal_quintic_hermite_aos must have three elements: {y, y', y''}");
 418         }
 419         if (dx <= 0)
 420         {
 421             throw std::domain_error("dx > 0 is required.");
 422         }
 423
 424         for (auto & datum : data_)
 425         {
 426             datum[1] *= dx;
 427             datum[2] *= (dx*dx/2);
 428         }
 429     }
 430
 431
 432     inline Real operator()(Real x) const
 433     {
 434         const Real xf = x0_ + (data_.size()-1)/inv_dx_;
 435         if  (x < x0_ || x > xf)
 436         {
 437             std::ostringstream oss;
 438             oss.precision(std::numeric_limits<Real>::digits10+3);
 439             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
 440                 << x0_ << ", " << xf << "]";
 441             throw std::domain_error(oss.str());
 442         }
 443         if (x == xf)
 444         {
 445             return data_.back()[0];
 446         }
 447         return this->unchecked_evaluation(x);
 448     }
 449
 450     inline Real unchecked_evaluation(Real x) const
 451     {
 452         using std::floor;
 453         Real s = (x-x0_)*inv_dx_;
 454         Real ii = floor(s);
 455         auto i = static_cast<decltype(data_.size())>(ii);
 456         Real t = s - ii;
 457         if (t == 0)
 458         {
 459             return data_[i][0];
 460         }
 461
 462         Real y0 = data_[i][0];
 463         Real dy0 = data_[i][1];
 464         Real d2y0 = data_[i][2];
 465         Real y1 = data_[i+1][0];
 466         Real dy1 = data_[i+1][1];
 467         Real d2y1 = data_[i+1][2];
 468
 469         Real y = (1 - t*t*t*(10 + t*(-15 + 6*t)))*y0;
 470         y += t*(1 + t*t*(-6 + t*(8 - 3*t)))*dy0;
 471         y += t*t*(1 + t*(-3 + t*(3 - t)))*d2y0;
 472         y += t*t*t*((1 + t*(-2 + t))*d2y1 + (-4 + t*(7 - 3*t))*dy1 + (10 + t*(-15 + 6*t))*y1);
 473         return y;
 474     }
 475
 476     inline Real prime(Real x) const
 477     {
 478         const Real xf = x0_ + (data_.size()-1)/inv_dx_;
 479         if  (x < x0_ || x > xf)
 480         {
 481             std::ostringstream oss;
 482             oss.precision(std::numeric_limits<Real>::digits10+3);
 483             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
 484                 << x0_ << ", " << xf << "]";
 485             throw std::domain_error(oss.str());
 486         }
 487         if (x == xf)
 488         {
 489             return data_.back()[1]*inv_dx_;
 490         }
 491
 492         return this->unchecked_prime(x);
 493     }
 494
 495     inline Real unchecked_prime(Real x) const
 496     {
 497         using std::floor;
 498         Real s = (x-x0_)*inv_dx_;
 499         Real ii = floor(s);
 500         auto i = static_cast<decltype(data_.size())>(ii);
 501         Real t = s - ii;
 502         if (t == 0)
 503         {
 504             return data_[i][1]*inv_dx_;
 505         }
 506
 507
 508         Real y0 = data_[i][0];
 509         Real y1 = data_[i+1][0];
 510         Real v0 = data_[i][1];
 511         Real v1 = data_[i+1][1];
 512         Real a0 = data_[i][2];
 513         Real a1 = data_[i+1][2];
 514
 515         Real dy = 30*t*t*(1 - 2*t + t*t)*(y1-y0);
 516         dy += (1-18*t*t + 32*t*t*t - 15*t*t*t*t)*v0 - t*t*(12 - 28*t + 15*t*t)*v1;
 517         dy += t*((2 - 9*t + 12*t*t - 5*t*t*t)*a0 + t*(3 - 8*t + 5*t*t)*a1);
 518         return dy*inv_dx_;
 519     }
 520
 521     inline Real double_prime(Real x) const
 522     {
 523         const Real xf = x0_ + (data_.size()-1)/inv_dx_;
 524         if  (x < x0_ || x > xf)
 525         {
 526             std::ostringstream oss;
 527             oss.precision(std::numeric_limits<Real>::digits10+3);
 528             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
 529                 << x0_ << ", " << xf << "]";
 530             throw std::domain_error(oss.str());
 531         }
 532         if (x == xf)
 533         {
 534             return data_.back()[2]*2*inv_dx_*inv_dx_;
 535         }
 536
 537         return this->unchecked_double_prime(x);
 538     }
 539
 540     inline Real unchecked_double_prime(Real x) const
 541     {
 542         using std::floor;
 543         Real s = (x-x0_)*inv_dx_;
 544         Real ii = floor(s);
 545         auto i = static_cast<decltype(data_.size())>(ii);
 546         Real t = s - ii;
 547         if (t == 0) {
 548             return data_[i][2]*2*inv_dx_*inv_dx_;
 549         }
 550         Real y0 = data_[i][0];
 551         Real dy0 = data_[i][1];
 552         Real d2y0 = data_[i][2];
 553         Real y1 = data_[i+1][0];
 554         Real dy1 = data_[i+1][1];
 555         Real d2y1 = data_[i+1][2];
 556
 557         Real d2ydx2 = 60*t*(1 - 3*t + 2*t*t)*(y1 - y0)*inv_dx_*inv_dx_;
 558         d2ydx2 += (12*t)*((-3 + 8*t - 5*t*t)*dy0 - (2 - 7*t + 5*t*t)*dy1);
 559         d2ydx2 += (1 - 9*t + 18*t*t - 10*t*t*t)*d2y0*(2*inv_dx_*inv_dx_) + t*(3 - 12*t + 10*t*t)*d2y1*(2*inv_dx_*inv_dx_);
 560         return d2ydx2;
 561     }
 562
 563     int64_t bytes() const
 564     {
 565         return data_.size()*data_[0].size()*sizeof(Real) + 2*sizeof(Real);
 566     }
 567
 568     std::pair<Real, Real> domain() const
 569     {
 570         Real xf = x0_ + (data_.size()-1)/inv_dx_;
 571         return {x0_, xf};
 572     }
 573
 574 private:
 575     RandomAccessContainer data_;
 576     Real x0_;
 577     Real inv_dx_;
 578 };
 579
 580 }
 581 }
 582 }
 583 }
 584 #endif