# HG changeset patch
# User Alexander Barth
# Date 1204786718 18000
# Node ID 1507a9d6df401ad8907f75ab0637d219f8e9a523
# Parent  b1368dc9e719f9e9748fc5ccc09dbb66b75f749f
__lin_interpn__.cc: handle decreasing coordinate values

diff -r b1368dc9e719 -r 1507a9d6df40 src/ChangeLog
--- a/src/ChangeLog	Wed Mar 05 04:58:54 2008 -0500
+++ b/src/ChangeLog	Thu Mar 06 01:58:38 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-02-27  John W. Eaton  <jwe@octave.org>
 
 	* oct-stream.cc (do_read): Stop reading if seek fails.
diff -r b1368dc9e719 -r 1507a9d6df40 src/DLD-FUNCTIONS/__lin_interpn__.cc
--- a/src/DLD-FUNCTIONS/__lin_interpn__.cc	Wed Mar 05 04:58:54 2008 -0500
+++ b/src/DLD-FUNCTIONS/__lin_interpn__.cc	Thu Mar 06 01:58:38 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