ceph/src/boost/boost/math/interpolators/detail/septic_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_SEPTIC_HERMITE_DETAIL_HPP
   8 #define BOOST_MATH_INTERPOLATORS_DETAIL_SEPTIC_HERMITE_DETAIL_HPP
   9 #include <algorithm>
  10 #include <stdexcept>
  11 #include <sstream>
  12 #include <cmath>
  13
  14 namespace boost::math::interpolators::detail {
  15
  16 template<class RandomAccessContainer>
  17 class septic_hermite_detail {
  18 public:
  19     using Real = typename RandomAccessContainer::value_type;
  20     septic_hermite_detail(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2, RandomAccessContainer && d3ydx3)
  21     : x_{std::move(x)}, y_{std::move(y)}, dydx_{std::move(dydx)}, d2ydx2_{std::move(d2ydx2)}, d3ydx3_{std::move(d3ydx3)}
  22     {
  23         if (x_.size() != y_.size())
  24         {
  25             throw std::domain_error("Number of abscissas must = number of ordinates.");
  26         }
  27         if (x_.size() != dydx_.size())
  28         {
  29             throw std::domain_error("Numbers of derivatives must = number of abscissas.");
  30         }
  31         if (x_.size() != d2ydx2_.size())
  32         {
  33             throw std::domain_error("Number of second derivatives must equal number of abscissas.");
  34         }
  35         if (x_.size() != d3ydx3_.size())
  36         {
  37             throw std::domain_error("Number of third derivatives must equal number of abscissas.");
  38         }
  39
  40         if (x_.size() < 2)
  41         {
  42             throw std::domain_error("At least 2 abscissas are required.");
  43         }
  44         Real x0 = x_[0];
  45         for (decltype(x_.size()) i = 1; i < x_.size(); ++i)
  46         {
  47             Real x1 = x_[i];
  48             if (x1 <= x0)
  49             {
  50                 throw std::domain_error("Abscissas must be sorted in strictly increasing order x0 < x1 < ... < x_{n-1}");
  51             }
  52             x0 = x1;
  53         }
  54     }
  55
  56     void push_back(Real x, Real y, Real dydx, Real d2ydx2, Real d3ydx3)
  57     {
  58         using std::abs;
  59         using std::isnan;
  60         if (x <= x_.back()) {
  61              throw std::domain_error("Calling push_back must preserve the monotonicity of the x's");
  62         }
  63         x_.push_back(x);
  64         y_.push_back(y);
  65         dydx_.push_back(dydx);
  66         d2ydx2_.push_back(d2ydx2);
  67         d3ydx3_.push_back(d3ydx3);
  68     }
  69
  70     Real operator()(Real x) const
  71     {
  72         if  (x < x_[0] || x > x_.back())
  73         {
  74             std::ostringstream oss;
  75             oss.precision(std::numeric_limits<Real>::digits10+3);
  76             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
  77                 << x_[0] << ", " << x_.back() << "]";
  78             throw std::domain_error(oss.str());
  79         }
  80         // t \in [0, 1)
  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 dx = (x1-x0);
  91         Real t = (x-x0)/dx;
  92
  93         // See:
  94         // http://seisweb.usask.ca/classes/GEOL481/2017/Labs/interpolation_utilities_matlab/shermite.m
  95         Real t2 = t*t;
  96         Real t3 = t2*t;
  97         Real t4 = t3*t;
  98         Real dx2 = dx*dx/2;
  99         Real dx3 = dx2*dx/3;
 100
 101         Real s = t4*(-35 + t*(84 + t*(-70 + 20*t)));
 102         Real z4 = -s;
 103         Real z0 = s + 1;
 104         Real z1 = t*(1 + t3*(-20 + t*(45 + t*(-36 + 10*t))));
 105         Real z2 = t2*(1 + t2*(-10 + t*(20 + t*(-15 + 4*t))));
 106         Real z3 = t3*(1 + t*(-4 + t*(6 + t*(-4 + t))));
 107         Real z5 = t4*(-15 + t*(39 + t*(-34 + 10*t)));
 108         Real z6 = t4*(5 + t*(-14 + t*(13 - 4*t)));
 109         Real z7 = t4*(-1 + t*(3 + t*(-3+t)));
 110
 111         Real y0 = y_[i];
 112         Real y1 = y_[i+1];
 113         // Velocity:
 114         Real v0 = dydx_[i];
 115         Real v1 = dydx_[i+1];
 116         // Acceleration:
 117         Real a0 = d2ydx2_[i];
 118         Real a1 = d2ydx2_[i+1];
 119         // Jerk:
 120         Real j0 = d3ydx3_[i];
 121         Real j1 = d3ydx3_[i+1];
 122
 123         return z0*y0 + z4*y1 + (z1*v0 + z5*v1)*dx + (z2*a0 + z6*a1)*dx2 + (z3*j0 + z7*j1)*dx3;
 124     }
 125
 126     Real prime(Real x) const
 127     {
 128         if  (x < x_[0] || x > x_.back())
 129         {
 130             std::ostringstream oss;
 131             oss.precision(std::numeric_limits<Real>::digits10+3);
 132             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
 133                 << x_[0] << ", " << x_.back() << "]";
 134             throw std::domain_error(oss.str());
 135         }
 136         if (x == x_.back())
 137         {
 138             return dydx_.back();
 139         }
 140
 141         auto it = std::upper_bound(x_.begin(), x_.end(), x);
 142         auto i = std::distance(x_.begin(), it) -1;
 143         Real x0 = *(it-1);
 144         Real x1 = *it;
 145         Real y0 = y_[i];
 146         Real y1 = y_[i+1];
 147         Real v0 = dydx_[i];
 148         Real v1 = dydx_[i+1];
 149         Real a0 = d2ydx2_[i];
 150         Real a1 = d2ydx2_[i+1];
 151         Real j0 = d3ydx3_[i];
 152         Real j1 = d3ydx3_[i+1];
 153         Real dx = x1 - x0;
 154         Real t = (x-x0)/dx;
 155         Real t2 = t*t;
 156         Real t3 = t2*t;
 157         Real z0 = 140*t3*(1 + t*(-3 + t*(3 - t)));
 158         Real z1 = 1 + t3*(-80 + t*(225 + t*(-216 + 70*t)));
 159         Real z2 = t3*(-60 + t*(195 + t*(-204 + 70*t)));
 160         Real z3 = 1 + t2*(-20 + t*(50 + t*(-45 + 14*t)));
 161         Real z4 = t2*(10 + t*(-35 + t*(39 - 14*t)));
 162         Real z5 = 3 + t*(-16 + t*(30 + t*(-24 + 7*t)));
 163         Real z6 = t*(-4 + t*(15 + t*(-18 + 7*t)));
 164
 165         Real dydx = z0*(y1-y0)/dx;
 166         dydx += z1*v0 + z2*v1;
 167         dydx += (x-x0)*(z3*a0 + z4*a1);
 168         dydx += (x-x0)*(x-x0)*(z5*j0 + z6*j1)/6;
 169         return dydx;
 170     }
 171
 172     inline Real double_prime(Real x) const
 173     {
 174         return std::numeric_limits<Real>::quiet_NaN();
 175     }
 176
 177     friend std::ostream& operator<<(std::ostream & os, const septic_hermite_detail & m)
 178     {
 179         os << "(x,y,y') = {";
 180         for (size_t i = 0; i < m.x_.size() - 1; ++i) {
 181             os << "(" << m.x_[i] << ", " << m.y_[i] << ", " << m.dydx_[i] << ", " << m.d2ydx2_[i] <<  ", " << m.d3ydx3_[i] << "),  ";
 182         }
 183         auto n = m.x_.size()-1;
 184         os << "(" << m.x_[n] << ", " << m.y_[n] << ", " << m.dydx_[n] << ", " << m.d2ydx2_[n] << m.d3ydx3_[n] << ")}";
 185         return os;
 186     }
 187
 188     int64_t bytes()
 189     {
 190         return 5*x_.size()*sizeof(Real) + 5*sizeof(x_);
 191     }
 192
 193     std::pair<Real, Real> domain() const
 194     {
 195         return {x_.front(), x_.back()};
 196     }
 197
 198 private:
 199     RandomAccessContainer x_;
 200     RandomAccessContainer y_;
 201     RandomAccessContainer dydx_;
 202     RandomAccessContainer d2ydx2_;
 203     RandomAccessContainer d3ydx3_;
 204 };
 205
 206 template<class RandomAccessContainer>
 207 class cardinal_septic_hermite_detail {
 208 public:
 209     using Real = typename RandomAccessContainer::value_type;
 210     cardinal_septic_hermite_detail(RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2, RandomAccessContainer && d3ydx3, Real x0, Real dx)
 211     : y_{std::move(y)}, dy_{std::move(dydx)}, d2y_{std::move(d2ydx2)}, d3y_{std::move(d3ydx3)}, x0_{x0}, inv_dx_{1/dx}
 212     {
 213         if (y_.size() != dy_.size())
 214         {
 215             throw std::domain_error("Numbers of derivatives must = number of ordinates.");
 216         }
 217         if (y_.size() != d2y_.size())
 218         {
 219             throw std::domain_error("Number of second derivatives must equal number of ordinates.");
 220         }
 221         if (y_.size() != d3y_.size())
 222         {
 223             throw std::domain_error("Number of third derivatives must equal number of ordinates.");
 224         }
 225         if (y_.size() < 2)
 226         {
 227             throw std::domain_error("At least 2 abscissas are required.");
 228         }
 229
 230         if (dx <= 0)
 231         {
 232             throw std::domain_error("dx > 0 is required.");
 233         }
 234
 235         for (auto & dy : dy_)
 236         {
 237             dy *= dx;
 238         }
 239         for (auto & d2y : d2y_)
 240         {
 241             d2y *= (dx*dx/2);
 242         }
 243         for (auto & d3y : d3y_)
 244         {
 245             d3y *= (dx*dx*dx/6);
 246         }
 247
 248     }
 249
 250     inline Real operator()(Real x) const
 251     {
 252         Real xf = x0_ + (y_.size()-1)/inv_dx_;
 253         if  (x < x0_ || x > xf)
 254         {
 255             std::ostringstream oss;
 256             oss.precision(std::numeric_limits<Real>::digits10+3);
 257             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
 258                 << x0_ << ", " << xf << "]";
 259             throw std::domain_error(oss.str());
 260         }
 261         if (x == xf)
 262         {
 263             return y_.back();
 264         }
 265         return this->unchecked_evaluation(x);
 266     }
 267
 268     inline Real unchecked_evaluation(Real x) const {
 269         using std::floor;
 270         Real s3 = (x-x0_)*inv_dx_;
 271         Real ii = floor(s3);
 272         auto i = static_cast<decltype(y_.size())>(ii);
 273         Real t = s3 - ii;
 274         if (t == 0) {
 275             return y_[i];
 276         }
 277         // See:
 278         // http://seisweb.usask.ca/classes/GEOL481/2017/Labs/interpolation_utilities_matlab/shermite.m
 279         Real t2 = t*t;
 280         Real t3 = t2*t;
 281         Real t4 = t3*t;
 282
 283         Real s = t4*(-35 + t*(84 + t*(-70 + 20*t)));
 284         Real z4 = -s;
 285         Real z0 = s + 1;
 286         Real z1 = t*(1 + t3*(-20 + t*(45 + t*(-36+10*t))));
 287         Real z2 = t2*(1 + t2*(-10 + t*(20 + t*(-15+4*t))));
 288         Real z3 = t3*(1 + t*(-4+t*(6+t*(-4+t))));
 289         Real z5 = t4*(-15 + t*(39 + t*(-34 + 10*t)));
 290         Real z6 = t4*(5 + t*(-14 + t*(13-4*t)));
 291         Real z7 = t4*(-1 + t*(3+t*(-3+t)));
 292
 293         Real y0 = y_[i];
 294         Real y1 = y_[i+1];
 295         Real dy0 = dy_[i];
 296         Real dy1 = dy_[i+1];
 297         Real a0 = d2y_[i];
 298         Real a1 = d2y_[i+1];
 299         Real j0 = d3y_[i];
 300         Real j1 = d3y_[i+1];
 301
 302         return z0*y0 + z1*dy0 + z2*a0 + z3*j0 + z4*y1 + z5*dy1 + z6*a1 + z7*j1;
 303     }
 304
 305     inline Real prime(Real x) const
 306     {
 307         Real xf = x0_ + (y_.size()-1)/inv_dx_;
 308         if  (x < x0_ || x > xf)
 309         {
 310             std::ostringstream oss;
 311             oss.precision(std::numeric_limits<Real>::digits10+3);
 312             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
 313                 << x0_ << ", " << xf << "]";
 314             throw std::domain_error(oss.str());
 315         }
 316         if (x == xf)
 317         {
 318             return dy_.back()/inv_dx_;
 319         }
 320
 321         return this->unchecked_prime(x);
 322     }
 323
 324     inline Real unchecked_prime(Real x) const
 325     {
 326         using std::floor;
 327         Real s3 = (x-x0_)*inv_dx_;
 328         Real ii = floor(s3);
 329         auto i = static_cast<decltype(y_.size())>(ii);
 330         Real t = s3 - ii;
 331         if (t==0)
 332         {
 333             return dy_[i]/inv_dx_;
 334         }
 335
 336         Real y0 = y_[i];
 337         Real y1 = y_[i+1];
 338         Real dy0 = dy_[i];
 339         Real dy1 = dy_[i+1];
 340         Real a0 = d2y_[i];
 341         Real a1 = d2y_[i+1];
 342         Real j0 = d3y_[i];
 343         Real j1 = d3y_[i+1];
 344         Real t2 = t*t;
 345         Real t3 = t2*t;
 346         Real z0 = 140*t3*(1 + t*(-3 + t*(3 - t)));
 347         Real z1 = 1 + t3*(-80 + t*(225 + t*(-216 + 70*t)));
 348         Real z2 = t3*(-60 + t*(195 + t*(-204 + 70*t)));
 349         Real z3 = 1 + t2*(-20 + t*(50 + t*(-45 + 14*t)));
 350         Real z4 = t2*(10 + t*(-35 + t*(39 - 14*t)));
 351         Real z5 = 3 + t*(-16 + t*(30 + t*(-24 + 7*t)));
 352         Real z6 = t*(-4 + t*(15 + t*(-18 + 7*t)));
 353
 354         Real dydx = z0*(y1-y0)*inv_dx_;
 355         dydx += (z1*dy0 + z2*dy1)*inv_dx_;
 356         dydx += 2*t*(z3*a0 + z4*a1)*inv_dx_;
 357         dydx += t*t*(z5*j0 + z6*j1);
 358         return dydx;
 359     }
 360
 361     inline Real double_prime(Real x) const
 362     {
 363         Real xf = x0_ + (y_.size()-1)/inv_dx_;
 364         if  (x < x0_ || x > xf)
 365         {
 366             std::ostringstream oss;
 367             oss.precision(std::numeric_limits<Real>::digits10+3);
 368             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
 369                 << x0_ << ", " << xf << "]";
 370             throw std::domain_error(oss.str());
 371         }
 372         if (x == xf)
 373         {
 374             return d2y_.back()*2*inv_dx_*inv_dx_;
 375         }
 376
 377         return this->unchecked_double_prime(x);
 378     }
 379
 380     inline Real unchecked_double_prime(Real x) const
 381     {
 382         using std::floor;
 383         Real s3 = (x-x0_)*inv_dx_;
 384         Real ii = floor(s3);
 385         auto i = static_cast<decltype(y_.size())>(ii);
 386         Real t = s3 - ii;
 387         if (t==0)
 388         {
 389             return d2y_[i]*2*inv_dx_*inv_dx_;
 390         }
 391
 392         Real y0 = y_[i];
 393         Real y1 = y_[i+1];
 394         Real dy0 = dy_[i];
 395         Real dy1 = dy_[i+1];
 396         Real a0 = d2y_[i];
 397         Real a1 = d2y_[i+1];
 398         Real j0 = d3y_[i];
 399         Real j1 = d3y_[i+1];
 400         Real t2 = t*t;
 401
 402         Real z0 = 420*t2*(1 + t*(-4 + t*(5 - 2*t)));
 403         Real z1 = 60*t2*(-4 + t*(15 + t*(-18 + 7*t)));
 404         Real z2 = 60*t2*(-3 + t*(13 + t*(-17 + 7*t)));
 405         Real z3 = (1 + t2*(-60 + t*(200 + t*(-225 + 84*t))));
 406         Real z4 = t2*(30 + t*(-140 + t*(195 - 84*t)));
 407         Real z5 = t*(1 + t*(-8 + t*(20 + t*(-20 + 7*t))));
 408         Real z6 = t2*(-2 + t*(10 + t*(-15 + 7*t)));
 409
 410         Real d2ydx2 = z0*(y1-y0)*inv_dx_*inv_dx_;
 411         d2ydx2 += (z1*dy0 + z2*dy1)*inv_dx_*inv_dx_;
 412         d2ydx2 += (z3*a0 + z4*a1)*2*inv_dx_*inv_dx_;
 413         d2ydx2 += 6*(z5*j0 + z6*j1)/(inv_dx_*inv_dx_);
 414
 415         return d2ydx2;
 416     }
 417
 418     int64_t bytes() const
 419     {
 420         return 4*y_.size()*sizeof(Real) + 2*sizeof(Real) + 4*sizeof(y_);
 421     }
 422
 423     std::pair<Real, Real> domain() const
 424     {
 425         return {x0_, x0_ + (y_.size()-1)/inv_dx_};
 426     }
 427
 428 private:
 429     RandomAccessContainer y_;
 430     RandomAccessContainer dy_;
 431     RandomAccessContainer d2y_;
 432     RandomAccessContainer d3y_;
 433     Real x0_;
 434     Real inv_dx_;
 435 };
 436
 437
 438 template<class RandomAccessContainer>
 439 class cardinal_septic_hermite_detail_aos {
 440 public:
 441     using Point = typename RandomAccessContainer::value_type;
 442     using Real = typename Point::value_type;
 443     cardinal_septic_hermite_detail_aos(RandomAccessContainer && dat, Real x0, Real dx)
 444     : data_{std::move(dat)}, x0_{x0}, inv_dx_{1/dx}
 445     {
 446         if (data_.size() < 2) {
 447             throw std::domain_error("At least 2 abscissas are required.");
 448         }
 449         if (data_[0].size() != 4) {
 450             throw std::domain_error("There must be 4 data items per struct.");
 451         }
 452
 453         for (auto & datum : data_)
 454         {
 455             datum[1] *= dx;
 456             datum[2] *= (dx*dx/2);
 457             datum[3] *= (dx*dx*dx/6);
 458         }
 459     }
 460
 461     inline Real operator()(Real x) const
 462     {
 463         Real xf = x0_ + (data_.size()-1)/inv_dx_;
 464         if  (x < x0_ || x > xf)
 465         {
 466             std::ostringstream oss;
 467             oss.precision(std::numeric_limits<Real>::digits10+3);
 468             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
 469                 << x0_ << ", " << xf << "]";
 470             throw std::domain_error(oss.str());
 471         }
 472         if (x == xf)
 473         {
 474             return data_.back()[0];
 475         }
 476         return this->unchecked_evaluation(x);
 477     }
 478
 479     inline Real unchecked_evaluation(Real x) const
 480     {
 481         using std::floor;
 482         Real s3 = (x-x0_)*inv_dx_;
 483         Real ii = floor(s3);
 484         auto i = static_cast<decltype(data_.size())>(ii);
 485         Real t = s3 - ii;
 486         if (t==0)
 487         {
 488             return data_[i][0];
 489         }
 490         Real t2 = t*t;
 491         Real t3 = t2*t;
 492         Real t4 = t3*t;
 493
 494         Real s = t4*(-35 + t*(84 + t*(-70 + 20*t)));
 495         Real z4 = -s;
 496         Real z0 = s + 1;
 497         Real z1 = t*(1 + t3*(-20 + t*(45 + t*(-36+10*t))));
 498         Real z2 = t2*(1 + t2*(-10 + t*(20 + t*(-15+4*t))));
 499         Real z3 = t3*(1 + t*(-4+t*(6+t*(-4+t))));
 500         Real z5 = t4*(-15 + t*(39 + t*(-34 + 10*t)));
 501         Real z6 = t4*(5 + t*(-14 + t*(13-4*t)));
 502         Real z7 = t4*(-1 + t*(3+t*(-3+t)));
 503
 504         Real y0 = data_[i][0];
 505         Real dy0 = data_[i][1];
 506         Real a0 = data_[i][2];
 507         Real j0 = data_[i][3];
 508         Real y1 = data_[i+1][0];
 509         Real dy1 = data_[i+1][1];
 510         Real a1 = data_[i+1][2];
 511         Real j1 = data_[i+1][3];
 512
 513         return z0*y0 + z1*dy0 + z2*a0 + z3*j0 + z4*y1 + z5*dy1 + z6*a1 + z7*j1;
 514     }
 515
 516     inline Real prime(Real x) const
 517     {
 518         Real xf = x0_ + (data_.size()-1)/inv_dx_;
 519         if  (x < x0_ || x > xf)
 520         {
 521             std::ostringstream oss;
 522             oss.precision(std::numeric_limits<Real>::digits10+3);
 523             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
 524                 << x0_ << ", " << xf << "]";
 525             throw std::domain_error(oss.str());
 526         }
 527         if (x == xf)
 528         {
 529             return data_.back()[1]*inv_dx_;
 530         }
 531
 532         return this->unchecked_prime(x);
 533     }
 534
 535     inline Real unchecked_prime(Real x) const
 536     {
 537         using std::floor;
 538         Real s3 = (x-x0_)*inv_dx_;
 539         Real ii = floor(s3);
 540         auto i = static_cast<decltype(data_.size())>(ii);
 541         Real t = s3 - ii;
 542         if (t == 0)
 543         {
 544             return data_[i][1]*inv_dx_;
 545         }
 546
 547         Real y0 = data_[i][0];
 548         Real y1 = data_[i+1][0];
 549         Real dy0 = data_[i][1];
 550         Real dy1 = data_[i+1][1];
 551         Real a0 = data_[i][2];
 552         Real a1 = data_[i+1][2];
 553         Real j0 = data_[i][3];
 554         Real j1 = data_[i+1][3];
 555         Real t2 = t*t;
 556         Real t3 = t2*t;
 557         Real z0 = 140*t3*(1 + t*(-3 + t*(3 - t)));
 558         Real z1 = 1 + t3*(-80 + t*(225 + t*(-216 + 70*t)));
 559         Real z2 = t3*(-60 + t*(195 + t*(-204 + 70*t)));
 560         Real z3 = 1 + t2*(-20 + t*(50 + t*(-45 + 14*t)));
 561         Real z4 = t2*(10 + t*(-35 + t*(39 - 14*t)));
 562         Real z5 = 3 + t*(-16 + t*(30 + t*(-24 + 7*t)));
 563         Real z6 = t*(-4 + t*(15 + t*(-18 + 7*t)));
 564
 565         Real dydx = z0*(y1-y0)*inv_dx_;
 566         dydx += (z1*dy0 + z2*dy1)*inv_dx_;
 567         dydx += 2*t*(z3*a0 + z4*a1)*inv_dx_;
 568         dydx += t*t*(z5*j0 + z6*j1);
 569         return dydx;
 570     }
 571
 572     inline Real double_prime(Real x) const
 573     {
 574         Real xf = x0_ + (data_.size()-1)/inv_dx_;
 575         if  (x < x0_ || x > xf)
 576         {
 577             std::ostringstream oss;
 578             oss.precision(std::numeric_limits<Real>::digits10+3);
 579             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
 580                 << x0_ << ", " << xf << "]";
 581             throw std::domain_error(oss.str());
 582         }
 583         if (x == xf)
 584         {
 585             return data_.back()[2]*2*inv_dx_*inv_dx_;
 586         }
 587
 588         return this->unchecked_double_prime(x);
 589     }
 590
 591     inline Real unchecked_double_prime(Real x) const
 592     {
 593         using std::floor;
 594         Real s3 = (x-x0_)*inv_dx_;
 595         Real ii = floor(s3);
 596         auto i = static_cast<decltype(data_.size())>(ii);
 597         Real t = s3 - ii;
 598         if (t == 0)
 599         {
 600             return data_[i][2]*2*inv_dx_*inv_dx_;
 601         }
 602         Real y0 = data_[i][0];
 603         Real y1 = data_[i+1][0];
 604         Real dy0 = data_[i][1];
 605         Real dy1 = data_[i+1][1];
 606         Real a0 = data_[i][2];
 607         Real a1 = data_[i+1][2];
 608         Real j0 = data_[i][3];
 609         Real j1 = data_[i+1][3];
 610         Real t2 = t*t;
 611
 612         Real z0 = 420*t2*(1 + t*(-4 + t*(5 - 2*t)));
 613         Real z1 = 60*t2*(-4 + t*(15 + t*(-18 + 7*t)));
 614         Real z2 = 60*t2*(-3 + t*(13 + t*(-17 + 7*t)));
 615         Real z3 = (1 + t2*(-60 + t*(200 + t*(-225 + 84*t))));
 616         Real z4 = t2*(30 + t*(-140 + t*(195 - 84*t)));
 617         Real z5 = t*(1 + t*(-8 + t*(20 + t*(-20 + 7*t))));
 618         Real z6 = t2*(-2 + t*(10 + t*(-15 + 7*t)));
 619
 620         Real d2ydx2 = z0*(y1-y0)*inv_dx_*inv_dx_;
 621         d2ydx2 += (z1*dy0 + z2*dy1)*inv_dx_*inv_dx_;
 622         d2ydx2 += (z3*a0 + z4*a1)*2*inv_dx_*inv_dx_;
 623         d2ydx2 += 6*(z5*j0 + z6*j1)/(inv_dx_*inv_dx_);
 624
 625         return d2ydx2;
 626     }
 627
 628     int64_t bytes() const
 629     {
 630         return data_.size()*data_[0].size()*sizeof(Real) + 2*sizeof(Real) + sizeof(data_);
 631     }
 632
 633     std::pair<Real, Real> domain() const
 634     {
 635         return {x0_, x0_ + (data_.size() -1)/inv_dx_};
 636     }
 637
 638 private:
 639     RandomAccessContainer data_;
 640     Real x0_;
 641     Real inv_dx_;
 642 };
 643
 644
 645 }
 646 #endif