Warning, /include/Geant4/tools/spline is written in an unsupported language. File is not indexed.
0001 // Copyright (C) 2010, Guy Barrand. All rights reserved.
0002 // See the file tools.license for terms.
0003
0004 #ifndef tools_spline
0005 #define tools_spline
0006
0007 // From Federico Carminati code found in root-6.08.06/TSpline.h, TSpline.cxx.
0008
0009 #include "mnmx"
0010 #include <cstddef>
0011 #include <vector>
0012 #include <ostream>
0013 #include <cmath>
0014
0015 namespace tools {
0016 namespace spline {
0017
0018 class base_poly {
0019 public:
0020 base_poly():fX(0),fY(0) {}
0021 base_poly(double x,double y):fX(x),fY(y) {}
0022 virtual ~base_poly(){}
0023 public:
0024 base_poly(base_poly const &a_from):fX(a_from.fX),fY(a_from.fY) {}
0025 base_poly& operator=(base_poly const &a_from) {
0026 if(this==&a_from) return *this;
0027 fX = a_from.fX;
0028 fY = a_from.fY;
0029 return *this;
0030 }
0031 public:
0032 const double& X() const {return fX;}
0033 const double& Y() const {return fY;}
0034 double &X() {return fX;}
0035 double &Y() {return fY;}
0036 protected:
0037 double fX; // abscissa
0038 double fY; // constant term
0039 };
0040
0041 class cubic_poly : public base_poly {
0042 public:
0043 cubic_poly():fB(0), fC(0), fD(0) {}
0044 cubic_poly(double x, double y, double b, double c, double d):base_poly(x,y), fB(b), fC(c), fD(d) {}
0045 public:
0046 cubic_poly(cubic_poly const &a_from)
0047 :base_poly(a_from), fB(a_from.fB), fC(a_from.fC), fD(a_from.fD) {}
0048 cubic_poly& operator=(cubic_poly const &a_from) {
0049 if(this==&a_from) return *this;
0050 base_poly::operator=(a_from);
0051 fB = a_from.fB;
0052 fC = a_from.fC;
0053 fD = a_from.fD;
0054 return *this;
0055 }
0056 public:
0057 double &B() {return fB;}
0058 double &C() {return fC;}
0059 double &D() {return fD;}
0060 double eval(double x) const {double dx=x-fX;return (fY+dx*(fB+dx*(fC+dx*fD)));}
0061 protected:
0062 double fB; // first order expansion coefficient : fB*1! is the first derivative at x
0063 double fC; // second order expansion coefficient : fC*2! is the second derivative at x
0064 double fD; // third order expansion coefficient : fD*3! is the third derivative at x
0065 };
0066
0067 class quintic_poly : public base_poly {
0068 public:
0069 quintic_poly():fB(0), fC(0), fD(0), fE(0), fF(0) {}
0070 quintic_poly(double x, double y, double b, double c, double d, double e, double f)
0071 :base_poly(x,y), fB(b), fC(c), fD(d), fE(e), fF(f) {}
0072 public:
0073 quintic_poly(quintic_poly const &a_from)
0074 :base_poly(a_from)
0075 ,fB(a_from.fB),fC(a_from.fC),fD(a_from.fD),fE(a_from.fE),fF(a_from.fF) {}
0076 quintic_poly& operator=(quintic_poly const &a_from) {
0077 if(this==&a_from) return *this;
0078 base_poly::operator=(a_from);
0079 fB = a_from.fB;
0080 fC = a_from.fC;
0081 fD = a_from.fD;
0082 fE = a_from.fE;
0083 fF = a_from.fF;
0084 return *this;
0085 }
0086 public:
0087 double &B() {return fB;}
0088 double &C() {return fC;}
0089 double &D() {return fD;}
0090 double &E() {return fE;}
0091 double &F() {return fF;}
0092 double eval(double x) const {double dx=x-fX;return (fY+dx*(fB+dx*(fC+dx*(fD+dx*(fE+dx*fF)))));}
0093 protected:
0094 double fB; // first order expansion coefficient : fB*1! is the first derivative at x
0095 double fC; // second order expansion coefficient : fC*2! is the second derivative at x
0096 double fD; // third order expansion coefficient : fD*3! is the third derivative at x
0097 double fE; // fourth order expansion coefficient : fE*4! is the fourth derivative at x
0098 double fF; // fifth order expansion coefficient : fF*5! is the fifth derivative at x
0099 };
0100
0101 ////////////////////////////////////////////////////////////////////////////////////
0102 ////////////////////////////////////////////////////////////////////////////////////
0103 ////////////////////////////////////////////////////////////////////////////////////
0104 class base_spline {
0105 protected:
0106 base_spline(std::ostream& a_out):m_out(a_out), fDelta(-1), fXmin(0), fXmax(0), fNp(0), fKstep(false) {}
0107 public:
0108 base_spline(std::ostream& a_out,double delta, double xmin, double xmax, size_t np, bool step)
0109 :m_out(a_out),fDelta(delta), fXmin(xmin),fXmax(xmax), fNp(np), fKstep(step)
0110 {}
0111 virtual ~base_spline() {}
0112 protected:
0113 base_spline(const base_spline& a_from)
0114 :m_out(a_from.m_out)
0115 ,fDelta(a_from.fDelta),fXmin(a_from.fXmin),fXmax(a_from.fXmax),fNp(a_from.fNp),fKstep(a_from.fKstep) {}
0116 base_spline& operator=(const base_spline& a_from) {
0117 if(this==&a_from) return *this;
0118 fDelta=a_from.fDelta;
0119 fXmin=a_from.fXmin;
0120 fXmax=a_from.fXmax;
0121 fNp=a_from.fNp;
0122 fKstep=a_from.fKstep;
0123 return *this;
0124 }
0125 protected:
0126 std::ostream& m_out;
0127 double fDelta; // Distance between equidistant knots
0128 double fXmin; // Minimum value of abscissa
0129 double fXmax; // Maximum value of abscissa
0130 size_t fNp; // Number of knots
0131 bool fKstep; // True of equidistant knots
0132 };
0133
0134
0135 //////////////////////////////////////////////////////////////////////////
0136 // //
0137 // cubic //
0138 // //
0139 // Class to create third splines to interpolate knots //
0140 // Arbitrary conditions can be introduced for first and second //
0141 // derivatives at beginning and ending points //
0142 // //
0143 //////////////////////////////////////////////////////////////////////////
0144
0145 class cubic : public base_spline {
0146 protected:
0147 cubic(std::ostream& a_out) : base_spline(a_out) , fPoly(0), fValBeg(0), fValEnd(0), fBegCond(-1), fEndCond(-1) {}
0148 public:
0149 cubic(std::ostream& a_out,size_t a_n,double a_x[], double a_y[], double a_valbeg = 0, double a_valend = 0)
0150 :base_spline(a_out,-1,0,0,a_n,false)
0151 ,fValBeg(a_valbeg), fValEnd(a_valend), fBegCond(0), fEndCond(0)
0152 {
0153 if(!a_n) {
0154 m_out << "tools::spline::cubic : a_np is null." << std::endl;
0155 return;
0156 }
0157 fXmin = a_x[0];
0158 fXmax = a_x[a_n-1];
0159 fPoly.resize(a_n);
0160 for (size_t i=0; i<a_n; ++i) {
0161 fPoly[i].X() = a_x[i];
0162 fPoly[i].Y() = a_y[i];
0163 }
0164 build_coeff(); // Build the spline coefficients
0165 }
0166 public:
0167 cubic(const cubic& a_from)
0168 :base_spline(a_from)
0169 ,fPoly(a_from.fPoly),fValBeg(a_from.fValBeg),fValEnd(a_from.fValEnd),fBegCond(a_from.fBegCond),fEndCond(a_from.fEndCond)
0170 {}
0171 cubic& operator=(const cubic& a_from) {
0172 if(this==&a_from) return *this;
0173 base_spline::operator=(a_from);
0174 fPoly = a_from.fPoly;
0175 fValBeg=a_from.fValBeg;
0176 fValEnd=a_from.fValEnd;
0177 fBegCond=a_from.fBegCond;
0178 fEndCond=a_from.fEndCond;
0179 return *this;
0180 }
0181 public:
0182 double eval(double x) const {
0183 if(!fNp) return 0;
0184 // Eval this spline at x
0185 size_t klow = find_x(x);
0186 if ( (fNp > 1) && (klow >= (fNp-1))) klow = fNp-2; //see: https://savannah.cern.ch/bugs/?71651
0187 return fPoly[klow].eval(x);
0188 }
0189 protected:
0190 template<typename T>
0191 static int TMath_Nint(T x) {
0192 // Round to nearest integer. Rounds half integers to the nearest even integer.
0193 int i;
0194 if (x >= 0) {
0195 i = int(x + 0.5);
0196 if ( i & 1 && x + 0.5 == T(i) ) i--;
0197 } else {
0198 i = int(x - 0.5);
0199 if ( i & 1 && x - 0.5 == T(i) ) i++;
0200 }
0201 return i;
0202 }
0203 static int TMath_FloorNint(double x) { return TMath_Nint(::floor(x)); }
0204
0205 size_t find_x(double x) const {
0206 int klow=0, khig=int(fNp-1);
0207 //
0208 // If out of boundaries, extrapolate
0209 // It may be badly wrong
0210 if(x<=fXmin) klow=0;
0211 else if(x>=fXmax) klow=khig;
0212 else {
0213 if(fKstep) { // Equidistant knots, use histogramming :
0214 klow = TMath_FloorNint((x-fXmin)/fDelta);
0215 // Correction for rounding errors
0216 if (x < fPoly[klow].X())
0217 klow = mx<int>(klow-1,0);
0218 else if (klow < khig) {
0219 if (x > fPoly[klow+1].X()) ++klow;
0220 }
0221 } else {
0222 int khalf;
0223 //
0224 // Non equidistant knots, binary search
0225 while((khig-klow)>1) {
0226 khalf = (klow+khig)/2;
0227 if(x>fPoly[khalf].X()) klow=khalf;
0228 else khig=khalf;
0229 }
0230 //
0231 // This could be removed, sanity check
0232 if( (x<fPoly[klow].X()) || (fPoly[klow+1].X()<x) ) {
0233 m_out << "tools::spline::cubic::find_x : Binary search failed"
0234 << " x(" << klow << ") = " << fPoly[klow].X() << " < x= " << x
0235 << " < x(" << klow+1 << ") = " << fPoly[klow+1].X() << "."
0236 << "." << std::endl;
0237 }
0238 }
0239 }
0240 return klow;
0241 }
0242
0243 void build_coeff() {
0244 /// subroutine cubspl ( tau, c, n, ibcbeg, ibcend )
0245 /// from * a practical guide to splines * by c. de boor
0246 /// ************************ input ***************************
0247 /// n = number of data points. assumed to be .ge. 2.
0248 /// (tau(i), c(1,i), i=1,...,n) = abscissae and ordinates of the
0249 /// data points. tau is assumed to be strictly increasing.
0250 /// ibcbeg, ibcend = boundary condition indicators, and
0251 /// c(2,1), c(2,n) = boundary condition information. specifically,
0252 /// ibcbeg = 0 means no boundary condition at tau(1) is given.
0253 /// in this case, the not-a-knot condition is used, i.e. the
0254 /// jump in the third derivative across tau(2) is forced to
0255 /// zero, thus the first and the second cubic polynomial pieces
0256 /// are made to coincide.)
0257 /// ibcbeg = 1 means that the slope at tau(1) is made to equal
0258 /// c(2,1), supplied by input.
0259 /// ibcbeg = 2 means that the second derivative at tau(1) is
0260 /// made to equal c(2,1), supplied by input.
0261 /// ibcend = 0, 1, or 2 has analogous meaning concerning the
0262 /// boundary condition at tau(n), with the additional infor-
0263 /// mation taken from c(2,n).
0264 /// *********************** output **************************
0265 /// c(j,i), j=1,...,4; i=1,...,l (= n-1) = the polynomial coefficients
0266 /// of the cubic interpolating spline with interior knots (or
0267 /// joints) tau(2), ..., tau(n-1). precisely, in the interval
0268 /// (tau(i), tau(i+1)), the spline f is given by
0269 /// f(x) = c(1,i)+h*(c(2,i)+h*(c(3,i)+h*c(4,i)/3.)/2.)
0270 /// where h = x - tau(i). the function program *ppvalu* may be
0271 /// used to evaluate f or its derivatives from tau,c, l = n-1,
0272 /// and k=4.
0273
0274 int j, l;
0275 double divdf1,divdf3,dtau,g=0;
0276 // ***** a tridiagonal linear system for the unknown slopes s(i) of
0277 // f at tau(i), i=1,...,n, is generated and then solved by gauss elim-
0278 // ination, with s(i) ending up in c(2,i), all i.
0279 // c(3,.) and c(4,.) are used initially for temporary storage.
0280 l = int(fNp-1);
0281 // compute first differences of x sequence and store in C also,
0282 // compute first divided difference of data and store in D.
0283 {for (size_t m=1; m<fNp ; ++m) {
0284 fPoly[m].C() = fPoly[m].X() - fPoly[m-1].X();
0285 fPoly[m].D() = (fPoly[m].Y() - fPoly[m-1].Y())/fPoly[m].C();
0286 }}
0287 // construct first equation from the boundary condition, of the form
0288 // D[0]*s[0] + C[0]*s[1] = B[0]
0289 if(fBegCond==0) {
0290 if(fNp == 2) {
0291 // no condition at left end and n = 2.
0292 fPoly[0].D() = 1.;
0293 fPoly[0].C() = 1.;
0294 fPoly[0].B() = 2.*fPoly[1].D();
0295 } else {
0296 // not-a-knot condition at left end and n .gt. 2.
0297 fPoly[0].D() = fPoly[2].C();
0298 fPoly[0].C() = fPoly[1].C() + fPoly[2].C();
0299 fPoly[0].B() = ((fPoly[1].C()+2.*fPoly[0].C())*fPoly[1].D()*fPoly[2].C()+
0300 fPoly[1].C()*fPoly[1].C()*fPoly[2].D())/fPoly[0].C();
0301 }
0302 } else if (fBegCond==1) {
0303 // slope prescribed at left end.
0304 fPoly[0].B() = fValBeg;
0305 fPoly[0].D() = 1.;
0306 fPoly[0].C() = 0.;
0307 } else if (fBegCond==2) {
0308 // second derivative prescribed at left end.
0309 fPoly[0].D() = 2.;
0310 fPoly[0].C() = 1.;
0311 fPoly[0].B() = 3.*fPoly[1].D() - fPoly[1].C()/2.*fValBeg;
0312 }
0313 bool forward_gauss_elimination = true;
0314 if(fNp > 2) {
0315 // if there are interior knots, generate the corresp. equations and car-
0316 // ry out the forward pass of gauss elimination, after which the m-th
0317 // equation reads D[m]*s[m] + C[m]*s[m+1] = B[m].
0318 {for (int m=1; m<l; ++m) {
0319 g = -fPoly[m+1].C()/fPoly[m-1].D();
0320 fPoly[m].B() = g*fPoly[m-1].B() + 3.*(fPoly[m].C()*fPoly[m+1].D()+fPoly[m+1].C()*fPoly[m].D());
0321 fPoly[m].D() = g*fPoly[m-1].C() + 2.*(fPoly[m].C() + fPoly[m+1].C());
0322 }}
0323 // construct last equation from the second boundary condition, of the form
0324 // (-g*D[n-2])*s[n-2] + D[n-1]*s[n-1] = B[n-1]
0325 // if slope is prescribed at right end, one can go directly to back-
0326 // substitution, since c array happens to be set up just right for it
0327 // at this point.
0328 if(fEndCond == 0) {
0329 if (fNp > 3 || fBegCond != 0) {
0330 // not-a-knot and n .ge. 3, and either n.gt.3 or also not-a-knot at
0331 // left end point.
0332 g = fPoly[fNp-2].C() + fPoly[fNp-1].C();
0333 fPoly[fNp-1].B() = ((fPoly[fNp-1].C()+2.*g)*fPoly[fNp-1].D()*fPoly[fNp-2].C()
0334 + fPoly[fNp-1].C()*fPoly[fNp-1].C()*(fPoly[fNp-2].Y()-fPoly[fNp-3].Y())/fPoly[fNp-2].C())/g;
0335 g = -g/fPoly[fNp-2].D();
0336 fPoly[fNp-1].D() = fPoly[fNp-2].C();
0337 } else {
0338 // either (n=3 and not-a-knot also at left) or (n=2 and not not-a-
0339 // knot at left end point).
0340 fPoly[fNp-1].B() = 2.*fPoly[fNp-1].D();
0341 fPoly[fNp-1].D() = 1.;
0342 g = -1./fPoly[fNp-2].D();
0343 }
0344 } else if (fEndCond == 1) {
0345 fPoly[fNp-1].B() = fValEnd;
0346 forward_gauss_elimination = false;
0347 } else if (fEndCond == 2) {
0348 // second derivative prescribed at right endpoint.
0349 fPoly[fNp-1].B() = 3.*fPoly[fNp-1].D() + fPoly[fNp-1].C()/2.*fValEnd;
0350 fPoly[fNp-1].D() = 2.;
0351 g = -1./fPoly[fNp-2].D();
0352 }
0353 } else {
0354 if(fEndCond == 0) {
0355 if (fBegCond > 0) {
0356 // either (n=3 and not-a-knot also at left) or (n=2 and not not-a-
0357 // knot at left end point).
0358 fPoly[fNp-1].B() = 2.*fPoly[fNp-1].D();
0359 fPoly[fNp-1].D() = 1.;
0360 g = -1./fPoly[fNp-2].D();
0361 } else {
0362 // not-a-knot at right endpoint and at left endpoint and n = 2.
0363 fPoly[fNp-1].B() = fPoly[fNp-1].D();
0364 forward_gauss_elimination = false;
0365 }
0366 } else if(fEndCond == 1) {
0367 fPoly[fNp-1].B() = fValEnd;
0368 forward_gauss_elimination = false;
0369 } else if(fEndCond == 2) {
0370 // second derivative prescribed at right endpoint.
0371 fPoly[fNp-1].B() = 3.*fPoly[fNp-1].D() + fPoly[fNp-1].C()/2.*fValEnd;
0372 fPoly[fNp-1].D() = 2.;
0373 g = -1./fPoly[fNp-2].D();
0374 }
0375 }
0376 // complete forward pass of gauss elimination.
0377 if(forward_gauss_elimination) {
0378 fPoly[fNp-1].D() = g*fPoly[fNp-2].C() + fPoly[fNp-1].D();
0379 fPoly[fNp-1].B() = (g*fPoly[fNp-2].B() + fPoly[fNp-1].B())/fPoly[fNp-1].D();
0380 }
0381 // carry out back substitution
0382 j = l-1;
0383 do {
0384 fPoly[j].B() = (fPoly[j].B() - fPoly[j].C()*fPoly[j+1].B())/fPoly[j].D();
0385 --j;
0386 } while (j>=0);
0387 // ****** generate cubic coefficients in each interval, i.e., the deriv.s
0388 // at its left endpoint, from value and slope at its endpoints.
0389 for (size_t i=1; i<fNp; ++i) {
0390 dtau = fPoly[i].C();
0391 divdf1 = (fPoly[i].Y() - fPoly[i-1].Y())/dtau;
0392 divdf3 = fPoly[i-1].B() + fPoly[i].B() - 2.*divdf1;
0393 fPoly[i-1].C() = (divdf1 - fPoly[i-1].B() - divdf3)/dtau;
0394 fPoly[i-1].D() = (divdf3/dtau)/dtau;
0395 }
0396 }
0397
0398 protected:
0399 std::vector<cubic_poly> fPoly; //[fNp] Array of polynomial terms
0400 double fValBeg; // Initial value of first or second derivative
0401 double fValEnd; // End value of first or second derivative
0402 int fBegCond; // 0=no beg cond, 1=first derivative, 2=second derivative
0403 int fEndCond; // 0=no end cond, 1=first derivative, 2=second derivative
0404 };
0405
0406 //////////////////////////////////////////////////////////////////////////
0407 // //
0408 // quintic //
0409 // //
0410 // Class to create quintic natural splines to interpolate knots //
0411 // Arbitrary conditions can be introduced for first and second //
0412 // derivatives using double knots (see build_coeff) for more on this. //
0413 // Double knots are automatically introduced at ending points //
0414 // //
0415 //////////////////////////////////////////////////////////////////////////
0416 class quintic : public base_spline {
0417 protected:
0418 quintic(std::ostream& a_out):base_spline(a_out),fPoly() {}
0419 public:
0420 quintic(std::ostream& a_out,size_t a_n ,double a_x[], double a_y[])
0421 :base_spline(a_out,-1,0,0,a_n,false) {
0422 if(!a_n) {
0423 m_out << "tools::spline::quintic : a_np is null." << std::endl;
0424 return;
0425 }
0426 fXmin = a_x[0];
0427 fXmax = a_x[a_n-1];
0428 fPoly.resize(fNp);
0429 for (size_t i=0; i<a_n; ++i) {
0430 fPoly[i].X() = a_x[i];
0431 fPoly[i].Y() = a_y[i];
0432 }
0433 build_coeff(); // Build the spline coefficients.
0434 }
0435 public:
0436 quintic(const quintic& a_from):base_spline(a_from),fPoly(a_from.fPoly) {}
0437 quintic& operator=(const quintic& a_from) {
0438 if(this==&a_from) return *this;
0439 base_spline::operator=(a_from);
0440 fPoly = a_from.fPoly;
0441 return *this;
0442 }
0443 public:
0444 double eval(double x) const {if(!fNp) return 0;size_t klow=find_x(x);return fPoly[klow].eval(x);}
0445
0446 protected:
0447 size_t find_x(double x) const {
0448 int klow=0;
0449
0450 // If out of boundaries, extrapolate
0451 // It may be badly wrong
0452 if(x<=fXmin) klow=0;
0453 else if(x>=fXmax) klow=int(fNp-1);
0454 else {
0455 if(fKstep) { // Equidistant knots, use histogramming :
0456 klow = mn<int>(int((x-fXmin)/fDelta),int(fNp-1));
0457 } else {
0458 int khig=int(fNp-1);
0459 int khalf;
0460 // Non equidistant knots, binary search
0461 while((khig-klow)>1) {
0462 khalf = (klow+khig)/2;
0463 if(x>fPoly[khalf].X()) klow=khalf;
0464 else khig=khalf;
0465 }
0466 }
0467
0468 // This could be removed, sanity check
0469 if( (x<fPoly[klow].X()) || (fPoly[klow+1].X()<x) ) {
0470 m_out << "tools::spline::quintic::find_x : Binary search failed"
0471 << " x(" << klow << ") = " << fPoly[klow].X() << " < x= " << x
0472 << " < x(" << klow+1<< ") = " << fPoly[klow+1].X() << "."
0473 << std::endl;
0474 }
0475 }
0476 return klow;
0477 }
0478
0479 void build_coeff() {
0480 ////////////////////////////////////////////////////////////////////////////////
0481 /// algorithm 600, collected algorithms from acm.
0482 /// algorithm appeared in acm-trans. math. software, vol.9, no. 2,
0483 /// jun., 1983, p. 258-259.
0484 ///
0485 /// quintic computes the coefficients of a quintic natural quintic spli
0486 /// s(x) with knots x(i) interpolating there to given function values:
0487 /// s(x(i)) = y(i) for i = 1,2, ..., n.
0488 /// in each interval (x(i),x(i+1)) the spline function s(xx) is a
0489 /// polynomial of fifth degree:
0490 /// s(xx) = ((((f(i)*p+e(i))*p+d(i))*p+c(i))*p+b(i))*p+y(i) (*)
0491 /// = ((((-f(i)*q+e(i+1))*q-d(i+1))*q+c(i+1))*q-b(i+1))*q+y(i+1)
0492 /// where p = xx - x(i) and q = x(i+1) - xx.
0493 /// (note the first subscript in the second expression.)
0494 /// the different polynomials are pieced together so that s(x) and
0495 /// its derivatives up to s"" are continuous.
0496 ///
0497 /// input:
0498 ///
0499 /// n number of data points, (at least three, i.e. n > 2)
0500 /// x(1:n) the strictly increasing or decreasing sequence of
0501 /// knots. the spacing must be such that the fifth power
0502 /// of x(i+1) - x(i) can be formed without overflow or
0503 /// underflow of exponents.
0504 /// y(1:n) the prescribed function values at the knots.
0505 ///
0506 /// output:
0507 ///
0508 /// b,c,d,e,f the computed spline coefficients as in (*).
0509 /// (1:n) specifically
0510 /// b(i) = s'(x(i)), c(i) = s"(x(i))/2, d(i) = s"'(x(i))/6,
0511 /// e(i) = s""(x(i))/24, f(i) = s""'(x(i))/120.
0512 /// f(n) is neither used nor altered. the five arrays
0513 /// b,c,d,e,f must always be distinct.
0514 ///
0515 /// option:
0516 ///
0517 /// it is possible to specify values for the first and second
0518 /// derivatives of the spline function at arbitrarily many knots.
0519 /// this is done by relaxing the requirement that the sequence of
0520 /// knots be strictly increasing or decreasing. specifically:
0521 ///
0522 /// if x(j) = x(j+1) then s(x(j)) = y(j) and s'(x(j)) = y(j+1),
0523 /// if x(j) = x(j+1) = x(j+2) then in addition s"(x(j)) = y(j+2).
0524 ///
0525 /// note that s""(x) is discontinuous at a double knot and, in
0526 /// addition, s"'(x) is discontinuous at a triple knot. the
0527 /// subroutine assigns y(i) to y(i+1) in these cases and also to
0528 /// y(i+2) at a triple knot. the representation (*) remains
0529 /// valid in each open interval (x(i),x(i+1)). at a double knot,
0530 /// x(j) = x(j+1), the output coefficients have the following values:
0531 /// y(j) = s(x(j)) = y(j+1)
0532 /// b(j) = s'(x(j)) = b(j+1)
0533 /// c(j) = s"(x(j))/2 = c(j+1)
0534 /// d(j) = s"'(x(j))/6 = d(j+1)
0535 /// e(j) = s""(x(j)-0)/24 e(j+1) = s""(x(j)+0)/24
0536 /// f(j) = s""'(x(j)-0)/120 f(j+1) = s""'(x(j)+0)/120
0537 /// at a triple knot, x(j) = x(j+1) = x(j+2), the output
0538 /// coefficients have the following values:
0539 /// y(j) = s(x(j)) = y(j+1) = y(j+2)
0540 /// b(j) = s'(x(j)) = b(j+1) = b(j+2)
0541 /// c(j) = s"(x(j))/2 = c(j+1) = c(j+2)
0542 /// d(j) = s"'((x(j)-0)/6 d(j+1) = 0 d(j+2) = s"'(x(j)+0)/6
0543 /// e(j) = s""(x(j)-0)/24 e(j+1) = 0 e(j+2) = s""(x(j)+0)/24
0544 /// f(j) = s""'(x(j)-0)/120 f(j+1) = 0 f(j+2) = s""'(x(j)+0)/120
0545
0546 size_t i, _m;
0547 double pqqr, p, q, r, _s, t, u, v,
0548 b1, p2, p3, q2, q3, r2, pq, pr, qr;
0549
0550 if (fNp <= 2) return;
0551
0552 // coefficients of a positive definite, pentadiagonal matrix,
0553 // stored in D, E, F from 1 to n-3.
0554 _m = fNp-2;
0555 q = fPoly[1].X()-fPoly[0].X();
0556 r = fPoly[2].X()-fPoly[1].X();
0557 q2 = q*q;
0558 r2 = r*r;
0559 qr = q+r;
0560 fPoly[0].D() = fPoly[0].E() = 0;
0561 if (q) fPoly[1].D() = q*6.*q2/(qr*qr);
0562 else fPoly[1].D() = 0;
0563
0564 if (_m > 1) {
0565 for (i = 1; i < _m; ++i) {
0566 p = q;
0567 q = r;
0568 r = fPoly[i+2].X()-fPoly[i+1].X();
0569 p2 = q2;
0570 q2 = r2;
0571 r2 = r*r;
0572 pq = qr;
0573 qr = q+r;
0574 if (q) {
0575 q3 = q2*q;
0576 pr = p*r;
0577 pqqr = pq*qr;
0578 fPoly[i+1].D() = q3*6./(qr*qr);
0579 fPoly[i].D() += (q+q)*(pr*15.*pr+(p+r)*q
0580 *(pr* 20.+q2*7.)+q2*
0581 ((p2+r2)*8.+pr*21.+q2+q2))/(pqqr*pqqr);
0582 fPoly[i-1].D() += q3*6./(pq*pq);
0583 fPoly[i].E() = q2*(p*qr+pq*3.*(qr+r+r))/(pqqr*qr);
0584 fPoly[i-1].E() += q2*(r*pq+qr*3.*(pq+p+p))/(pqqr*pq);
0585 fPoly[i-1].F() = q3/pqqr;
0586 } else
0587 fPoly[i+1].D() = fPoly[i].E() = fPoly[i-1].F() = 0;
0588 }
0589 }
0590 if (r) fPoly[_m-1].D() += r*6.*r2/(qr*qr);
0591
0592 // First and second order divided differences of the given function
0593 // values, stored in b from 2 to n and in c from 3 to n
0594 // respectively. care is taken of double and triple knots.
0595 for (i = 1; i < fNp; ++i) {
0596 if (fPoly[i].X() != fPoly[i-1].X()) {
0597 fPoly[i].B() =
0598 (fPoly[i].Y()-fPoly[i-1].Y())/(fPoly[i].X()-fPoly[i-1].X());
0599 } else {
0600 fPoly[i].B() = fPoly[i].Y();
0601 fPoly[i].Y() = fPoly[i-1].Y();
0602 }
0603 }
0604 for (i = 2; i < fNp; ++i) {
0605 if (fPoly[i].X() != fPoly[i-2].X()) {
0606 fPoly[i].C() =
0607 (fPoly[i].B()-fPoly[i-1].B())/(fPoly[i].X()-fPoly[i-2].X());
0608 } else {
0609 fPoly[i].C() = fPoly[i].B()*.5;
0610 fPoly[i].B() = fPoly[i-1].B();
0611 }
0612 }
0613
0614 // Solve the linear system with c(i+2) - c(i+1) as right-hand side.
0615 if (_m > 1) {
0616 p=fPoly[0].C()=fPoly[_m-1].E()=fPoly[0].F()
0617 =fPoly[_m-2].F()=fPoly[_m-1].F()=0;
0618 fPoly[1].C() = fPoly[3].C()-fPoly[2].C();
0619 fPoly[1].D() = 1./fPoly[1].D();
0620
0621 if (_m > 2) {
0622 for (i = 2; i < _m; ++i) {
0623 q = fPoly[i-1].D()*fPoly[i-1].E();
0624 fPoly[i].D() = 1./(fPoly[i].D()-p*fPoly[i-2].F()-q*fPoly[i-1].E());
0625 fPoly[i].E() -= q*fPoly[i-1].F();
0626 fPoly[i].C() = fPoly[i+2].C()-fPoly[i+1].C()-p*fPoly[i-2].C()
0627 -q*fPoly[i-1].C();
0628 p = fPoly[i-1].D()*fPoly[i-1].F();
0629 }
0630 }
0631 }
0632
0633 fPoly[fNp-2].C() = fPoly[fNp-1].C() = 0;
0634 if (fNp > 3)
0635 for (i=fNp-3; i > 0; --i)
0636 fPoly[i].C() = (fPoly[i].C()-fPoly[i].E()*fPoly[i+1].C()
0637 -fPoly[i].F()*fPoly[i+2].C())*fPoly[i].D();
0638
0639 // Integrate the third derivative of s(x)
0640 _m = fNp-1;
0641 q = fPoly[1].X()-fPoly[0].X();
0642 r = fPoly[2].X()-fPoly[1].X();
0643 b1 = fPoly[1].B();
0644 q3 = q*q*q;
0645 qr = q+r;
0646 if (qr) {
0647 v = fPoly[1].C()/qr;
0648 t = v;
0649 } else
0650 v = t = 0;
0651 if (q) fPoly[0].F() = v/q;
0652 else fPoly[0].F() = 0;
0653 for (i = 1; i < _m; ++i) {
0654 p = q;
0655 q = r;
0656 if (i != _m-1) r = fPoly[i+2].X()-fPoly[i+1].X();
0657 else r = 0;
0658 p3 = q3;
0659 q3 = q*q*q;
0660 pq = qr;
0661 qr = q+r;
0662 _s = t;
0663 if (qr) t = (fPoly[i+1].C()-fPoly[i].C())/qr;
0664 else t = 0;
0665 u = v;
0666 v = t-_s;
0667 if (pq) {
0668 fPoly[i].F() = fPoly[i-1].F();
0669 if (q) fPoly[i].F() = v/q;
0670 fPoly[i].E() = _s*5.;
0671 fPoly[i].D() = (fPoly[i].C()-q*_s)*10;
0672 fPoly[i].C() =
0673 fPoly[i].D()*(p-q)+(fPoly[i+1].B()-fPoly[i].B()+(u-fPoly[i].E())*
0674 p3-(v+fPoly[i].E())*q3)/pq;
0675 fPoly[i].B() = (p*(fPoly[i+1].B()-v*q3)+q*(fPoly[i].B()-u*p3))/pq-p
0676 *q*(fPoly[i].D()+fPoly[i].E()*(q-p));
0677 } else {
0678 fPoly[i].C() = fPoly[i-1].C();
0679 fPoly[i].D() = fPoly[i].E() = fPoly[i].F() = 0;
0680 }
0681 }
0682
0683 // End points x(1) and x(n)
0684 p = fPoly[1].X()-fPoly[0].X();
0685 _s = fPoly[0].F()*p*p*p;
0686 fPoly[0].E() = fPoly[0].D() = 0;
0687 fPoly[0].C() = fPoly[1].C()-_s*10;
0688 fPoly[0].B() = b1-(fPoly[0].C()+_s)*p;
0689
0690 q = fPoly[fNp-1].X()-fPoly[fNp-2].X();
0691 t = fPoly[fNp-2].F()*q*q*q;
0692 fPoly[fNp-1].E() = fPoly[fNp-1].D() = 0;
0693 fPoly[fNp-1].C() = fPoly[fNp-2].C()+t*10;
0694 fPoly[fNp-1].B() += (fPoly[fNp-1].C()-t)*q;
0695 }
0696 protected:
0697 std::vector<quintic_poly> fPoly; //[fNp] Array of polynomial terms
0698 };
0699
0700 }}
0701
0702 #endif
