[Bug-gsl] Interpolation bug?

[Andrew W. Steiner]

Hello all,

On my machine, after running the code below, which attempts
to use gsl_interp_accel_find() to search data which looks like:

index value
0 +0.000000e+00
1 +1.000000e+00
2 +2.000000e+00
3 +3.000000e+00
4 +4.000000e+00
5 +5.000000e+00
6 +6.000000e+00
7 +7.000000e+00
8 +8.000000e+00
9 +9.000000e+00

The function accel_find(0.9) gives 0, and accel_find(1.0) also gives 0,
then accel_find(1.1) gives 1, and then accel_find(1.0) gives 1, instead of
zero, because the cache value is different. I think the second call to
accel_find() is in error, and should actually give 1 rather than 0. Of
course, this is floating point arithmetic, but possibly there is a modification
to accel_find() which will make more sense most of the time.

Consider changing line 64 of interpolation/accel.c from

  else if (x > xa[x_index + 1])

to

  else if (x >= xa[x_index + 1])

which gives the result I expect.

HTH,
Andrew

--------------------------------------------------------------------------

for(i=0;i<10;i++) {
  x[i]=((double)(i));
  cout << i << " " << x[i] << endl;
}
cout << endl;

gsl_interp_accel *a=gsl_interp_accel_alloc();
cout << gsl_interp_accel_find(a,x,10,0.9) << endl;
cout << gsl_interp_accel_find(a,x,10,1.0) << endl;
cout << gsl_interp_accel_find(a,x,10,1.1) << endl;
cout << gsl_interp_accel_find(a,x,10,1.0) << endl;