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);
}
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
