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