[Bug-gsl] bug report regarding gsl_diff_*

[Zegenhagen Stefan]

Dear all,

I think I've found a bug that makes the numerical differentiation
routines quite unusable in some cases.

As visible from the source file diff/diff.c, in all the differentiation
routines, you calculate a higher order derivative first to choose a
stepsize that minimizes the absolute error when calculating f'. As a
stepsize to calculate the higher order derivative, you always choose the
absolute value of

GSL_SQRT_DBL_EPSILON.

This method fails if the value x at which to calculate the derivative is
so large that (x + i*GSL_SQRT_DBL_EPSILON) == x holds due to the
restricted numerical precision. Then, some of the intermediate values in
the calculation will become infinite (because of division by zero) and
the final result for h is NaN in the calculation of f'. This is
reproducable.

I also think that the same holds if x is so small that the equation (x +
i*GSL_SQRT_DBL_EPSILON) == GSL_SQRT_DBL_EPSILON holds. Both conditions
may be fulfilled even for double precision values.


The solution could be to either a stepsize relative to x in the
calculation of the higher order derivative(something like 

        double h = x * GSL_SQRT_DBL_EPSILON;

in the function header). This solution should avoid the restricted
precision problems.


Sincerely yours,
-- 

----------------------------------------------------------------------
Dipl. Phys. Stefan Zegenhagen
Max-Planck-Institut fuer Plasmaphysik
- Teilinstitut Greifswald -
Wendelsteinstr. 1
D-17491 Greifswald, Germany
Tel. +49 (0)3834 88 1277   Fax:  +49 (0)3834 88 2509
e-mail: address@hidden
-----------------------------------------------------------------------