Bug in gsl_interp_accel_find()

[Thomas Weber <thomas.weber.mail@.........>]

Hi, 

sorry, just realized that there's a separate list for bug reports
(originally posted to
http://lists.gnu.org/archive/html/help-gsl/2008-09/msg00003.html)


Consider the attached example code; the last assertion fails. Note that
idx2 == 0, but idx3 == 1.
The reason for that is that gsl_interp_accel_find() checks only for
strictly smaller/bigger values. 

(Untested) patch attached.

	Thomas
#include <gsl/gsl_vector.h>
#include <gsl/gsl_interp.h>

#include <stdio.h>
#include <assert.h>

int main()
{
  const double v[3] = {-0.1, 0.0, 0.1};

  size_t idx1, idx2, idx3;
  const double x1 = -0.05, x2 = 0;

  gsl_interp_accel * a = gsl_interp_accel_alloc();

  idx1 = gsl_interp_accel_find(a, v, 3, x1);
  idx2 = gsl_interp_accel_find(a, v, 3, x2);
  idx3 = gsl_interp_bsearch(v, x2, 0, 2);

  printf("idx1: %zi \t idx2: %zi \t idx3: %zi\n", idx1, idx2, idx3);

  assert( x2 < v[idx2+1]);
  
  return 0;
}
diff --git a/interpolation/gsl_interp.h b/interpolation/gsl_interp.h
index 7718be0..cfc60ab 100644
--- a/interpolation/gsl_interp.h
+++ b/interpolation/gsl_interp.h
@@ -205,7 +205,7 @@ gsl_interp_accel_find(gsl_interp_accel * a, const double xa[], size_t len, doubl
     a->miss_count++;
     a->cache = gsl_interp_bsearch(xa, x, 0, x_index);
   }
-  else if(x > xa[x_index + 1]) {
+  else if(x >= xa[x_index + 1]) {
     a->miss_count++;
     a->cache = gsl_interp_bsearch(xa, x, x_index, len-1);
   }