[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
## [Bug-gsl] Interpolation bug?

**From**: |
Andrew W. Steiner |

**Subject**: |
[Bug-gsl] Interpolation bug? |

**Date**: |
Thu, 25 Sep 2008 10:51:43 -0400 |

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;

**[Bug-gsl] Interpolation bug?**,
*Andrew W. Steiner* **<=**