[833] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
| 26 | // |
---|
[850] | 27 | // $Id: G4DataInterpolation.cc,v 1.10 2008/03/13 09:35:56 gcosmo Exp $ |
---|
| 28 | // GEANT4 tag $Name: HEAD $ |
---|
[833] | 29 | // |
---|
| 30 | #include "G4DataInterpolation.hh" |
---|
| 31 | |
---|
| 32 | ////////////////////////////////////////////////////////////////////////////// |
---|
| 33 | // |
---|
| 34 | // Constructor for initializing of fArgument, fFunction and fNumber |
---|
| 35 | // data members |
---|
| 36 | |
---|
| 37 | G4DataInterpolation::G4DataInterpolation( G4double pX[], |
---|
| 38 | G4double pY[], |
---|
| 39 | G4int number ) |
---|
| 40 | : fArgument(new G4double[number]), |
---|
| 41 | fFunction(new G4double[number]), |
---|
| 42 | fSecondDerivative(0), |
---|
| 43 | fNumber(number) |
---|
| 44 | { |
---|
| 45 | for(G4int i=0;i<fNumber;i++) |
---|
| 46 | { |
---|
| 47 | fArgument[i] = pX[i] ; |
---|
| 48 | fFunction[i] = pY[i] ; |
---|
| 49 | } |
---|
| 50 | } |
---|
| 51 | |
---|
| 52 | //////////////////////////////////////////////////////////////////////////// |
---|
| 53 | // |
---|
| 54 | // Constructor for cubic spline interpolation. It creates the array |
---|
| 55 | // fSecondDerivative[0,...fNumber-1] which is used in this interpolation by |
---|
| 56 | // the function |
---|
| 57 | |
---|
| 58 | |
---|
| 59 | G4DataInterpolation::G4DataInterpolation( G4double pX[], |
---|
| 60 | G4double pY[], |
---|
| 61 | G4int number, |
---|
| 62 | G4double pFirstDerStart, |
---|
| 63 | G4double pFirstDerFinish ) |
---|
| 64 | : fArgument(new G4double[number]), |
---|
| 65 | fFunction(new G4double[number]), |
---|
| 66 | fSecondDerivative(new G4double[number]), |
---|
| 67 | fNumber(number) |
---|
| 68 | { |
---|
| 69 | G4int i=0 ; |
---|
| 70 | G4double p=0.0, qn=0.0, sig=0.0, un=0.0 ; |
---|
| 71 | const G4double maxDerivative = 0.99e30 ; |
---|
| 72 | G4double* u = new G4double[fNumber - 1] ; |
---|
| 73 | |
---|
| 74 | for(i=0;i<fNumber;i++) |
---|
| 75 | { |
---|
| 76 | fArgument[i] = pX[i] ; |
---|
| 77 | fFunction[i] = pY[i] ; |
---|
| 78 | } |
---|
| 79 | if(pFirstDerStart > maxDerivative) |
---|
| 80 | { |
---|
| 81 | fSecondDerivative[0] = 0.0 ; |
---|
| 82 | u[0] = 0.0 ; |
---|
| 83 | } |
---|
| 84 | else |
---|
| 85 | { |
---|
| 86 | fSecondDerivative[0] = -0.5 ; |
---|
| 87 | u[0] = (3.0/(fArgument[1]-fArgument[0])) |
---|
| 88 | * ((fFunction[1]-fFunction[0])/(fArgument[1]-fArgument[0]) |
---|
| 89 | - pFirstDerStart) ; |
---|
| 90 | } |
---|
| 91 | |
---|
| 92 | // Decomposition loop for tridiagonal algorithm. fSecondDerivative[i] |
---|
| 93 | // and u[i] are used for temporary storage of the decomposed factors. |
---|
| 94 | |
---|
| 95 | for(i=1;i<fNumber-1;i++) |
---|
| 96 | { |
---|
| 97 | sig = (fArgument[i]-fArgument[i-1])/(fArgument[i+1]-fArgument[i-1]) ; |
---|
| 98 | p = sig*fSecondDerivative[i-1] + 2.0 ; |
---|
| 99 | fSecondDerivative[i] = (sig - 1.0)/p ; |
---|
| 100 | u[i] = (fFunction[i+1]-fFunction[i])/(fArgument[i+1]-fArgument[i]) - |
---|
| 101 | (fFunction[i]-fFunction[i-1])/(fArgument[i]-fArgument[i-1]) ; |
---|
| 102 | u[i] =(6.0*u[i]/(fArgument[i+1]-fArgument[i-1]) - sig*u[i-1])/p ; |
---|
| 103 | } |
---|
| 104 | if(pFirstDerFinish > maxDerivative) |
---|
| 105 | { |
---|
| 106 | qn = 0.0 ; |
---|
| 107 | un = 0.0 ; |
---|
| 108 | } |
---|
| 109 | else |
---|
| 110 | { |
---|
| 111 | qn = 0.5 ; |
---|
| 112 | un = (3.0/(fArgument[fNumber-1]-fArgument[fNumber-2])) |
---|
| 113 | * (pFirstDerFinish - (fFunction[fNumber-1]-fFunction[fNumber-2]) |
---|
| 114 | / (fArgument[fNumber-1]-fArgument[fNumber-2])) ; |
---|
| 115 | } |
---|
| 116 | fSecondDerivative[fNumber-1] = (un - qn*u[fNumber-2])/ |
---|
| 117 | (qn*fSecondDerivative[fNumber-2] + 1.0) ; |
---|
| 118 | |
---|
| 119 | // The backsubstitution loop for the triagonal algorithm of solving |
---|
| 120 | // a linear system of equations. |
---|
| 121 | |
---|
| 122 | for(G4int k=fNumber-2;k>=0;k--) |
---|
| 123 | { |
---|
| 124 | fSecondDerivative[k] = fSecondDerivative[k]*fSecondDerivative[k+1] + u[k]; |
---|
| 125 | } |
---|
| 126 | delete[] u ; |
---|
| 127 | } |
---|
| 128 | |
---|
| 129 | ///////////////////////////////////////////////////////////////////////////// |
---|
| 130 | // |
---|
| 131 | // Destructor deletes dynamically created arrays for data members: fArgument, |
---|
| 132 | // fFunction and fSecondDerivative, all have dimension of fNumber |
---|
| 133 | |
---|
| 134 | G4DataInterpolation::~G4DataInterpolation() |
---|
| 135 | { |
---|
| 136 | delete [] fArgument ; |
---|
| 137 | delete [] fFunction ; |
---|
| 138 | if(fSecondDerivative) { delete [] fSecondDerivative; } |
---|
| 139 | } |
---|
| 140 | |
---|
| 141 | ///////////////////////////////////////////////////////////////////////////// |
---|
| 142 | // |
---|
| 143 | // This function returns the value P(pX), where P(x) is polynom of fNumber-1 |
---|
| 144 | // degree such that P(fArgument[i]) = fFunction[i], for i = 0, ..., fNumber-1. |
---|
| 145 | // This is Lagrange's form of interpolation and it is based on Neville's |
---|
| 146 | // algorithm |
---|
| 147 | |
---|
| 148 | G4double |
---|
| 149 | G4DataInterpolation::PolynomInterpolation(G4double pX, |
---|
| 150 | G4double& deltaY ) const |
---|
| 151 | { |
---|
| 152 | G4int i=0, j=1, k=0 ; |
---|
| 153 | G4double mult=0.0, difi=0.0, deltaLow=0.0, deltaUp=0.0, cd=0.0, y=0.0 ; |
---|
| 154 | G4double* c = new G4double[fNumber] ; |
---|
| 155 | G4double* d = new G4double[fNumber] ; |
---|
| 156 | G4double diff = std::fabs(pX-fArgument[0]) ; |
---|
| 157 | for(i=0;i<fNumber;i++) |
---|
| 158 | { |
---|
| 159 | difi = std::fabs(pX-fArgument[i]) ; |
---|
| 160 | if(difi <diff) |
---|
| 161 | { |
---|
| 162 | k = i ; |
---|
| 163 | diff = difi ; |
---|
| 164 | } |
---|
| 165 | c[i] = fFunction[i] ; |
---|
| 166 | d[i] = fFunction[i] ; |
---|
| 167 | } |
---|
| 168 | y = fFunction[k--] ; |
---|
| 169 | for(j=1;j<fNumber;j++) |
---|
| 170 | { |
---|
| 171 | for(i=0;i<fNumber-j;i++) |
---|
| 172 | { |
---|
| 173 | deltaLow = fArgument[i] - pX ; |
---|
| 174 | deltaUp = fArgument[i+j] - pX ; |
---|
| 175 | cd = c[i+1] - d[i] ; |
---|
| 176 | mult = deltaLow - deltaUp ; |
---|
| 177 | if (!(mult != 0.0)) |
---|
| 178 | { |
---|
| 179 | G4Exception("G4DataInterpolation::PolynomInterpolation()", |
---|
| 180 | "Error", FatalException, "Coincident nodes !") ; |
---|
| 181 | } |
---|
| 182 | mult = cd/mult ; |
---|
| 183 | d[i] = deltaUp*mult ; |
---|
| 184 | c[i] = deltaLow*mult ; |
---|
| 185 | } |
---|
| 186 | y += (deltaY = (2*k < (fNumber - j -1) ? c[k+1] : d[k--] )) ; |
---|
| 187 | } |
---|
| 188 | delete[] c ; |
---|
| 189 | delete[] d ; |
---|
| 190 | |
---|
| 191 | return y ; |
---|
| 192 | } |
---|
| 193 | |
---|
| 194 | //////////////////////////////////////////////////////////////////////////// |
---|
| 195 | // |
---|
| 196 | // Given arrays fArgument[0,..,fNumber-1] and fFunction[0,..,fNumber-1], this |
---|
| 197 | // function calculates an array of coefficients. The coefficients don't provide |
---|
| 198 | // usually (fNumber>10) better accuracy for polynom interpolation, as compared |
---|
| 199 | // with PolynomInterpolation function. They could be used instead for derivate |
---|
| 200 | // calculations and some other applications. |
---|
| 201 | |
---|
| 202 | void |
---|
| 203 | G4DataInterpolation::PolIntCoefficient( G4double cof[]) const |
---|
| 204 | { |
---|
| 205 | G4int i=0, j=0 ; |
---|
| 206 | G4double factor=fNumber, reducedY=0.0, mult=1.0 ; |
---|
| 207 | G4double* tempArgument = new G4double[fNumber] ; |
---|
| 208 | |
---|
| 209 | for(i=0;i<fNumber;i++) |
---|
| 210 | { |
---|
| 211 | tempArgument[i] = cof[i] = 0.0 ; |
---|
| 212 | } |
---|
| 213 | tempArgument[fNumber-1] = -fArgument[0] ; |
---|
| 214 | |
---|
| 215 | for(i=1;i<fNumber;i++) |
---|
| 216 | { |
---|
| 217 | for(j=fNumber-1-i;j<fNumber-1;j++) |
---|
| 218 | { |
---|
| 219 | tempArgument[j] -= fArgument[i]*tempArgument[j+1] ; |
---|
| 220 | } |
---|
| 221 | tempArgument[fNumber-1] -= fArgument[i] ; |
---|
| 222 | } |
---|
| 223 | for(i=0;i<fNumber;i++) |
---|
| 224 | { |
---|
| 225 | factor = fNumber ; |
---|
| 226 | for(j=fNumber-1;j>=1;j--) |
---|
| 227 | { |
---|
| 228 | factor = j*tempArgument[j] + factor*fArgument[i] ; |
---|
| 229 | } |
---|
| 230 | reducedY = fFunction[i]/factor ; |
---|
| 231 | mult = 1.0 ; |
---|
| 232 | for(j=fNumber-1;j>=0;j--) |
---|
| 233 | { |
---|
| 234 | cof[j] += mult*reducedY ; |
---|
| 235 | mult = tempArgument[j] + mult*fArgument[i] ; |
---|
| 236 | } |
---|
| 237 | } |
---|
| 238 | delete[] tempArgument ; |
---|
| 239 | } |
---|
| 240 | |
---|
| 241 | ///////////////////////////////////////////////////////////////////////////// |
---|
| 242 | // |
---|
| 243 | // The function returns diagonal rational function (Bulirsch and Stoer |
---|
| 244 | // algorithm of Neville type) Pn(x)/Qm(x) where P and Q are polynoms. |
---|
| 245 | // Tests showed the method is not stable and hasn't advantage if compared |
---|
| 246 | // with polynomial interpolation ?! |
---|
| 247 | |
---|
| 248 | G4double |
---|
| 249 | G4DataInterpolation::RationalPolInterpolation(G4double pX, |
---|
| 250 | G4double& deltaY ) const |
---|
| 251 | { |
---|
| 252 | G4int i=0, j=1, k=0 ; |
---|
| 253 | const G4double tolerance = 1.6e-24 ; |
---|
| 254 | G4double mult=0.0, difi=0.0, cd=0.0, y=0.0, cof=0.0 ; |
---|
| 255 | G4double* c = new G4double[fNumber] ; |
---|
| 256 | G4double* d = new G4double[fNumber] ; |
---|
| 257 | G4double diff = std::fabs(pX-fArgument[0]) ; |
---|
| 258 | for(i=0;i<fNumber;i++) |
---|
| 259 | { |
---|
| 260 | difi = std::fabs(pX-fArgument[i]) ; |
---|
| 261 | if (!(difi != 0.0)) |
---|
| 262 | { |
---|
| 263 | y = fFunction[i] ; |
---|
| 264 | deltaY = 0.0 ; |
---|
| 265 | delete[] c ; |
---|
| 266 | delete[] d ; |
---|
| 267 | return y ; |
---|
| 268 | } |
---|
| 269 | else if(difi < diff) |
---|
| 270 | { |
---|
| 271 | k = i ; |
---|
| 272 | diff = difi ; |
---|
| 273 | } |
---|
| 274 | c[i] = fFunction[i] ; |
---|
| 275 | d[i] = fFunction[i] + tolerance ; // to prevent rare zero/zero cases |
---|
| 276 | } |
---|
| 277 | y = fFunction[k--] ; |
---|
| 278 | for(j=1;j<fNumber;j++) |
---|
| 279 | { |
---|
| 280 | for(i=0;i<fNumber-j;i++) |
---|
| 281 | { |
---|
| 282 | cd = c[i+1] - d[i] ; |
---|
| 283 | difi = fArgument[i+j] - pX ; |
---|
| 284 | cof = (fArgument[i] - pX)*d[i]/difi ; |
---|
| 285 | mult = cof - c[i+1] ; |
---|
| 286 | if (!(mult != 0.0)) // function to be interpolated has pole at pX |
---|
| 287 | { |
---|
| 288 | G4Exception("G4DataInterpolation::RationalPolInterpolation()", |
---|
| 289 | "Error", FatalException, "Coincident nodes !") ; |
---|
| 290 | } |
---|
| 291 | mult = cd/mult ; |
---|
| 292 | d[i] = c[i+1]*mult ; |
---|
| 293 | c[i] = cof*mult ; |
---|
| 294 | } |
---|
| 295 | y += (deltaY = (2*k < (fNumber - j - 1) ? c[k+1] : d[k--] )) ; |
---|
| 296 | } |
---|
| 297 | delete[] c ; |
---|
| 298 | delete[] d ; |
---|
| 299 | |
---|
| 300 | return y ; |
---|
| 301 | } |
---|
| 302 | |
---|
| 303 | ///////////////////////////////////////////////////////////////////////////// |
---|
| 304 | // |
---|
| 305 | // Cubic spline interpolation in point pX for function given by the table: |
---|
| 306 | // fArgument, fFunction. The constructor, which creates fSecondDerivative, |
---|
| 307 | // must be called before. The function works optimal, if sequential calls |
---|
| 308 | // are in random values of pX. |
---|
| 309 | |
---|
| 310 | G4double |
---|
| 311 | G4DataInterpolation::CubicSplineInterpolation(G4double pX) const |
---|
| 312 | { |
---|
| 313 | G4int kLow=0, kHigh=fNumber-1, k=0 ; |
---|
| 314 | |
---|
| 315 | // Searching in the table by means of bisection method. |
---|
| 316 | // fArgument must be monotonic, either increasing or decreasing |
---|
| 317 | |
---|
| 318 | while((kHigh - kLow) > 1) |
---|
| 319 | { |
---|
| 320 | k = (kHigh + kLow) >> 1 ; // compute midpoint 'bisection' |
---|
| 321 | if(fArgument[k] > pX) |
---|
| 322 | { |
---|
| 323 | kHigh = k ; |
---|
| 324 | } |
---|
| 325 | else |
---|
| 326 | { |
---|
| 327 | kLow = k ; |
---|
| 328 | } |
---|
| 329 | } // kLow and kHigh now bracket the input value of pX |
---|
| 330 | G4double deltaHL = fArgument[kHigh] - fArgument[kLow] ; |
---|
| 331 | if (!(deltaHL != 0.0)) |
---|
| 332 | { |
---|
| 333 | G4Exception("G4DataInterpolation::CubicSplineInterpolation()", |
---|
| 334 | "Error", FatalException, "Bad fArgument input !") ; |
---|
| 335 | } |
---|
| 336 | G4double a = (fArgument[kHigh] - pX)/deltaHL ; |
---|
| 337 | G4double b = (pX - fArgument[kLow])/deltaHL ; |
---|
| 338 | |
---|
| 339 | // Final evaluation of cubic spline polynomial for return |
---|
| 340 | |
---|
| 341 | return a*fFunction[kLow] + b*fFunction[kHigh] + |
---|
| 342 | ((a*a*a - a)*fSecondDerivative[kLow] + |
---|
| 343 | (b*b*b - b)*fSecondDerivative[kHigh])*deltaHL*deltaHL/6.0 ; |
---|
| 344 | } |
---|
| 345 | |
---|
| 346 | /////////////////////////////////////////////////////////////////////////// |
---|
| 347 | // |
---|
| 348 | // Return cubic spline interpolation in the point pX which is located between |
---|
| 349 | // fArgument[index] and fArgument[index+1]. It is usually called in sequence |
---|
| 350 | // of known from external analysis values of index. |
---|
| 351 | |
---|
| 352 | G4double |
---|
| 353 | G4DataInterpolation::FastCubicSpline(G4double pX, |
---|
| 354 | G4int index) const |
---|
| 355 | { |
---|
| 356 | G4double delta = fArgument[index+1] - fArgument[index] ; |
---|
| 357 | if (!(delta != 0.0)) |
---|
| 358 | { |
---|
| 359 | G4Exception("G4DataInterpolation::FastCubicSpline()", |
---|
| 360 | "Error", FatalException, "Bad fArgument input !") ; |
---|
| 361 | } |
---|
| 362 | G4double a = (fArgument[index+1] - pX)/delta ; |
---|
| 363 | G4double b = (pX - fArgument[index])/delta ; |
---|
| 364 | |
---|
| 365 | // Final evaluation of cubic spline polynomial for return |
---|
| 366 | |
---|
| 367 | return a*fFunction[index] + b*fFunction[index+1] + |
---|
| 368 | ((a*a*a - a)*fSecondDerivative[index] + |
---|
| 369 | (b*b*b - b)*fSecondDerivative[index+1])*delta*delta/6.0 ; |
---|
| 370 | } |
---|
| 371 | |
---|
| 372 | //////////////////////////////////////////////////////////////////////////// |
---|
| 373 | // |
---|
| 374 | // Given argument pX, returns index k, so that pX bracketed by fArgument[k] |
---|
| 375 | // and fArgument[k+1] |
---|
| 376 | |
---|
| 377 | G4int |
---|
| 378 | G4DataInterpolation::LocateArgument(G4double pX) const |
---|
| 379 | { |
---|
| 380 | G4int kLow=-1, kHigh=fNumber, k=0 ; |
---|
| 381 | G4bool ascend=(fArgument[fNumber-1] >= fArgument[0]) ; |
---|
| 382 | while((kHigh - kLow) > 1) |
---|
| 383 | { |
---|
| 384 | k = (kHigh + kLow) >> 1 ; // compute midpoint 'bisection' |
---|
| 385 | if( (pX >= fArgument[k]) == ascend) |
---|
| 386 | { |
---|
| 387 | kLow = k ; |
---|
| 388 | } |
---|
| 389 | else |
---|
| 390 | { |
---|
| 391 | kHigh = k ; |
---|
| 392 | } |
---|
| 393 | } |
---|
| 394 | if (!(pX != fArgument[0])) |
---|
| 395 | { |
---|
| 396 | return 1 ; |
---|
| 397 | } |
---|
| 398 | else if (!(pX != fArgument[fNumber-1])) |
---|
| 399 | { |
---|
| 400 | return fNumber - 2 ; |
---|
| 401 | } |
---|
| 402 | else return kLow ; |
---|
| 403 | } |
---|
| 404 | |
---|
| 405 | ///////////////////////////////////////////////////////////////////////////// |
---|
| 406 | // |
---|
| 407 | // Given a value pX, returns a value 'index' such that pX is between |
---|
| 408 | // fArgument[index] and fArgument[index+1]. fArgument MUST BE MONOTONIC, |
---|
| 409 | // either increasing or decreasing. If index = -1 or fNumber, this indicates |
---|
| 410 | // that pX is out of range. The value index on input is taken as the initial |
---|
| 411 | // approximation for index on output. |
---|
| 412 | |
---|
| 413 | void |
---|
| 414 | G4DataInterpolation::CorrelatedSearch( G4double pX, |
---|
| 415 | G4int& index ) const |
---|
| 416 | { |
---|
| 417 | G4int kHigh=0, k=0, Increment=0 ; |
---|
| 418 | // ascend = true for ascending order of table, false otherwise |
---|
| 419 | G4bool ascend = (fArgument[fNumber-1] >= fArgument[0]) ; |
---|
| 420 | if(index < 0 || index > fNumber-1) |
---|
| 421 | { |
---|
| 422 | index = -1 ; |
---|
| 423 | kHigh = fNumber ; |
---|
| 424 | } |
---|
| 425 | else |
---|
| 426 | { |
---|
| 427 | Increment = 1 ; // What value would be the best ? |
---|
| 428 | if((pX >= fArgument[index]) == ascend) |
---|
| 429 | { |
---|
| 430 | if(index == fNumber -1) |
---|
| 431 | { |
---|
| 432 | index = fNumber ; |
---|
| 433 | return ; |
---|
| 434 | } |
---|
| 435 | kHigh = index + 1 ; |
---|
| 436 | while((pX >= fArgument[kHigh]) == ascend) |
---|
| 437 | { |
---|
| 438 | index = kHigh ; |
---|
| 439 | Increment += Increment ; // double the Increment |
---|
| 440 | kHigh = index + Increment ; |
---|
| 441 | if(kHigh > (fNumber - 1)) |
---|
| 442 | { |
---|
| 443 | kHigh = fNumber ; |
---|
| 444 | break ; |
---|
| 445 | } |
---|
| 446 | } |
---|
| 447 | } |
---|
| 448 | else |
---|
| 449 | { |
---|
| 450 | if(index == 0) |
---|
| 451 | { |
---|
| 452 | index = -1 ; |
---|
| 453 | return ; |
---|
| 454 | } |
---|
| 455 | kHigh = index-- ; |
---|
| 456 | while((pX < fArgument[index]) == ascend) |
---|
| 457 | { |
---|
| 458 | kHigh = index ; |
---|
| 459 | Increment <<= 1 ; // double the Increment |
---|
| 460 | if(Increment >= kHigh) |
---|
| 461 | { |
---|
| 462 | index = -1 ; |
---|
| 463 | break ; |
---|
| 464 | } |
---|
| 465 | else |
---|
| 466 | { |
---|
| 467 | index = kHigh - Increment ; |
---|
| 468 | } |
---|
| 469 | } |
---|
| 470 | } // Value bracketed |
---|
| 471 | } |
---|
| 472 | // final bisection searching |
---|
| 473 | |
---|
| 474 | while((kHigh - index) != 1) |
---|
| 475 | { |
---|
| 476 | k = (kHigh + index) >> 1 ; |
---|
| 477 | if((pX >= fArgument[k]) == ascend) |
---|
| 478 | { |
---|
| 479 | index = k ; |
---|
| 480 | } |
---|
| 481 | else |
---|
| 482 | { |
---|
| 483 | kHigh = k ; |
---|
| 484 | } |
---|
| 485 | } |
---|
| 486 | if (!(pX != fArgument[fNumber-1])) |
---|
| 487 | { |
---|
| 488 | index = fNumber - 2 ; |
---|
| 489 | } |
---|
| 490 | if (!(pX != fArgument[0])) |
---|
| 491 | { |
---|
| 492 | index = 0 ; |
---|
| 493 | } |
---|
| 494 | return ; |
---|
| 495 | } |
---|
| 496 | |
---|
| 497 | // |
---|
| 498 | // |
---|
| 499 | //////////////////////////////////////////////////////////////////////////// |
---|