Mercurial > octave
diff src/DLD-FUNCTIONS/__lin_interpn__.cc @ 7561:a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
author | Alexander Barth |
---|---|
date | Thu, 06 Mar 2008 01:56:55 -0500 |
parents | 93c65f2a5668 |
children | 82be108cc558 72830070a17b |
line wrap: on
line diff
--- a/src/DLD-FUNCTIONS/__lin_interpn__.cc Wed Mar 05 14:24:33 2008 -0500 +++ b/src/DLD-FUNCTIONS/__lin_interpn__.cc Thu Mar 06 01:56:55 2008 -0500 @@ -43,38 +43,77 @@ octave_idx_type lookup (const double *x, octave_idx_type n, double y) { - octave_idx_type j, j0, j1; + octave_idx_type j; - if (y > x[n-1] || y < x[0]) - return -1; + if (x[0] < x[n-1]) + { + // increasing x + + if (y > x[n-1] || y < x[0]) + return -1; #ifdef EXHAUSTIF - for (j = 0; j < n - 1; j++) - { - if (x[j] <= y && y <= x[j+1]) - return j; - } + for (j = 0; j < n - 1; j++) + { + if (x[j] <= y && y <= x[j+1]) + return j; + } #else - j0 = 0; - j1 = n - 1; + octave_idx_type j0 = 0; + octave_idx_type j1 = n - 1; + + while (true) + { + j = (j0+j1)/2; + + if (y <= x[j+1]) + { + if (x[j] <= y) + return j; - while (true) + j1 = j; + } + + if (x[j] <= y) + j0 = j; + } +#endif + } + else { - j = (j0+j1)/2; + // decreasing x + // previous code with x -> -x and y -> -y - if (y <= x[j+1]) + if (y > x[0] || y < x[n-1]) + return -1; + +#ifdef EXHAUSTIF + for (j = 0; j < n - 1; j++) { - if (x[j] <= y) + if (x[j+1] <= y && y <= x[j]) return j; - - j1 = j; } +#else + octave_idx_type j0 = 0; + octave_idx_type j1 = n - 1; - if (x[j] <= y) - j0 = j; + while (true) + { + j = (j0+j1)/2; + + if (y >= x[j+1]) + { + if (x[j] >= y) + return j; + + j1 = j; + } + + if (x[j] >= y) + j0 = j; + } +#endif } - -#endif } // n-dimensional linear interpolation