With a small correction to the value of r (= M_PI / M_E), your modified test program (gsl-secant64.c) gives a correct result on Windows x64 using the version of secant.c now in the GSL repository:

# i   x_i-r      x_i - x_p status
  0 -7.509e-003  4.822e-002 -2
  1  0.000e+000  7.509e-003 -2
  2  0.000e+000  0.000e+000 0
# Converged
# f(x_i) =  2.2204e-019

Have you tested this revised version of secant.c?

Sorry for not participating in the discussion
in the middle of the week. Let us return to
the original issue. The key to reproduce it is
the operation system (linux and it was specified
in the original report). Concerning the version,
the last commit to secant.c was in 2007.

The difficulty with reproducing is related to difference
of the default FPU precision on Linux and on Windows
(80 and 64 bit accordingly).

My problem is not connected to the case when
the function evaluates exactly to zero.
It may happen that for some x_0
f(x_0) < 0 and f(x_0 + DBL_EPSILON) > 0
even for a linear function. So any test
for zero would not help.

If the function may be modified as (attached)
double test_f (double x, void *params)
{
   (void) params; // unused
   return -M_E*x + M_PI + 0.001*GSL_DBL_EPSILON;
}
then the problem may be reproduced with -mpc64
compile option, so I expect the same should happen
on Windows as well

$ cc -mpc64 -g -Wall -W   gsl-secant64.c -lgsl -lgslcblas -lm -o gsl-secant64
# i   x_i-r      x_i - x_p status
  0  2.830e-01  4.822e-02 -2
  1  2.905e-01  7.509e-03 -2
gsl: secant.c:104: ERROR: derivative value is not finite
Default GSL error handler invoked.

If the patch from the file file #28779 is applied
$ cc -mpc64 -g -Wall -W   gsl-secant64.c ../gsl-roots/libgsl-roots.a  -lgsl -lgslcblas -lm -o gsl-secant64
$ ./gsl-secant64
# i   x_i-r      x_i - x_p status
  0  2.830e-01  4.822e-02 -2
  1  2.905e-01  7.509e-03 -2
  2  2.905e-01  0.000e+00 0
# Converged
# f(x_i) =  2.2204e-19

Without -mpc the results are essentially the same
(# f(x_i) =  6.7273e-18 last case)

$ gcc --version
gcc (Ubuntu/Linaro 4.7.3-1ubuntu1) 4.7.3
Linux, i686 kernel,
Intel Core 2 Duo processor.

(file #29299)

I'm not objecting to a quick fix but I don't wan't to put a lot of effort into making a contribution to GSL if we can't also make improvements.  And I certainly consider that an API that in common situations forces twice as many function evaluations as are necessary to find a root is in need of improvement (in fact I really thought that I must be wrong about this and I still hope that I am).

I am sorry if this is not the right place to discuss such wider issues but I am new here and I wanted to understand the design as in my experience I have found it quite dangerous to 'fix' things that I don't understand.

My experience with having a compatibility define in MPIR (a fork of GMP) is quite the opposite to your own.  But, as you suggest, that discussion is for another place.

Not if we have to maintain the existing API's since we would need s new item in the state for the functiopn value.  And if you are saying that we have to maintain the existing API's for ever more, we cannot expect to improve GSL in areas where experience shows that API improvments are desirable.  I like the Python approach where they are willing to break legacy code every few years in order to accomodate improvements built on real experience (i.e GSL v2.0 would allow API changes).

One way round the legacy issue is to have two returns from the solvers SOLVER_SUCCESS and SOLVER_CONTINUE and to map these both onto GSL_SUCCESS for legacy code and onto GSL_SUCCESS and GSL_CONTINUE otherwise (do we have a define for detecting a user desire to use a legacy API?).

I agree that one could design a more crisp API (like you suggest) that'll better handle conditions like this, but you can't make the solver iteration routine return something other than GSL_SUCCESS without breaking legacy code.  Returning GSL_CONTINUE is nonzero and you'd break existing logic looking for failure from the iteration routine.

To stay within the current API/ABI, shortcircuiting on f == 0.0 and returning GSL_SUCCESS may be the only feasible change.  Anything else that's backwards compatible come to mind?

As for duplicated evalution, gsl_root_fsolver_root returns the root without requiring an additional function evaluation.   The root doesn't need to be exposed directly to the user through the iteration API as an accessor is available to get it at near-zero cost.  Again, one could imagine a different API but this is what's in-place now.  You can see this routine exercised in the examples and the source code in fsolver.c.

> So, in my view, the real question raised by this issue is what should the solver do if it is entered at a position where the previous iteration has already found the root and the user has not terminated the iteration?

This feels like a usage error on the user's part.  How one can have secant_iterate defend against having already found the root when "the solver found the root" is a problem- and tolerance-specific statement.  See below, however, ...

> To do nothing risks an infinite loop. But to do what is done now indicates failure for what is really success! Maybe a warning return value such as "root already found'?

Maybe, when the root is obtained identically (i.e. state->f == 0 on entry to secant_iterate in secant.c), have the method immediately return GSL_SUCCESS without attempting to re-evaluate the function.  I say identically because only in that case should all sane tolerancing schemes agree that the root has been found.

Presumably the user not checking for success will either (a) notice an infinite loop in their code or (b) notice that their solve procedure hits the maximum number of iterations if they control for that.

Thanks for digging into this,
Rhys

The section within secant.c in which the root position, the function value and the derivative value are updated is:

  x_new = x - (f / df);

  f_new = GSL_FN_FDF_EVAL_F(fdf, x_new) ;
  df_new = (f_new - f) / (x_new - x) ;

  *root = x_new ;

Although it involves an additional multiplication, I think this may be safer:

  x_new = x - (f / df);

  f_new = GSL_FN_FDF_EVAL_F(fdf, x_new) ;
  df_new = df * (1.0 - f_new / f) ;

  *root = x_new;

since it avoids the potential problem of x_new and x being the same because (f / df) is too small to change x when it is subtracted.  Also, since the algorithm shouldn't presumably be entered with f = 0, the division f_new / f should not cause a divide by zero problem.

My bad. Of course, the relative change of root approximation
matters, not the absolute one. The updated patch is attached.

(file #28783)