--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -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
+++ b/scripts/general/interpn.m
@@ -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
+++ b/src/ChangeLog
@@ -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
+++ b/src/DLD-FUNCTIONS/__lin_interpn__.cc
@@ -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