In several places in the diff.c code there are lines like the following:
h = sqrt (GSL_SQRT_DBL_EPSILON / (2.0 * a2)); ... *result = (GSL_FN_EVAL (f, x + h) - GSL_FN_EVAL (f, x)) / h; There is an assumption here that the difference between (x+h) and x is exactly identical to h, which is true only in infinite precision. In finite precision there can be a small error introduced because x+h may not be exactly representable in finite precision and so a nearby number is used (rounded or truncated). It is very easy to remove this error in the following way: h = sqrt (GSL_SQRT_DBL_EPSILON / (2.0 * a2)); ... h = (x+h)-x; // ensure difference between (x+h) and x is exactly h *result = (GSL_FN_EVAL (f, x + h) - GSL_FN_EVAL (f, x)) / h; The issue is explained in detail here - http://www.nrbook.com/a/bookcpdf/c5-7.pdf _______________________________________________ Bug-gsl mailing list [email protected] http://lists.gnu.org/mailman/listinfo/bug-gsl
