Mercurial > octave-nkf
changeset 7561:a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
| author | Alexander Barth |
|---|---|
| date | Thu, 06 Mar 2008 01:56:55 -0500 |
| parents | 0ef0f9802a37 |
| children | c827f5673321 |
| files | scripts/ChangeLog scripts/general/interpn.m src/ChangeLog src/DLD-FUNCTIONS/__lin_interpn__.cc |
| diffstat | 4 files changed, 73 insertions(+), 20 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog Wed Mar 05 14:24:33 2008 -0500 +++ b/scripts/ChangeLog Thu Mar 06 01:56:55 2008 -0500 @@ -1,3 +1,7 @@ +2008-03-06 John W. Eaton <jwe@octave.org> + + * general/interpn.m: New test. + 2008-03-05 Ben Abbott <bpabbott@mac.com> * polynomial/roots.m: Catch Infs and/or NaNs.
--- a/scripts/general/interpn.m Wed Mar 05 14:24:33 2008 -0500 +++ b/scripts/general/interpn.m Thu Mar 06 01:56:55 2008 -0500 @@ -250,3 +250,8 @@ %! assert (interpn(x,y,z,f,x,y,z), f) %! assert (interpn(x,y,z,f,x,y,z,'nearest'), f) %! assert (interpn(x,y,z,f,x,y,z,'spline'), f) + +%!test +%! [x,y,z] = ndgrid(0:2); +%! f = x.^2+y.^2+z.^2; +%! assert (interpn(x,y,-z,f,1.5,1.5,-1.5), 7.5)
--- a/src/ChangeLog Wed Mar 05 14:24:33 2008 -0500 +++ b/src/ChangeLog Thu Mar 06 01:56:55 2008 -0500 @@ -1,3 +1,8 @@ +2008-03-06 Alexander Barth <barth.alexander@gmail.com> + + * DLD-FUNCTIONS/__lin_interpn__.cc (lookup): + Handle decreasing coordinate values. + 2008-03-05 Jaroslav Hajek <highegg@gmail.com> * DLD-FUNCTIONS/chol.cc (Fcholupdate): Adjust code to meet
--- 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