diff --git a/roots/bisection.c b/roots/bisection.c
index 9665b736..6905b221 100644
--- a/roots/bisection.c
+++ b/roots/bisection.c
@@ -41,6 +41,7 @@ typedef struct
 bisection_state_t;
 
 static int bisection_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper);
+static int bisection_init2 (void * vstate, gsl_function * f, double * root, double x_lower, double f_lower, double x_upper, double f_upper);
 static int bisection_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper);
 
 static int
@@ -67,6 +68,25 @@ bisection_init (void * vstate, gsl_function * f, double * root, double x_lower,
 
 }
 
+static int
+bisection_init2 (void * vstate, gsl_function * f, double * root, double x_lower, double f_lower, double x_upper, double f_upper)
+{
+  bisection_state_t * state = (bisection_state_t *) vstate;
+
+  *root = 0.5 * (x_lower + x_upper) ;
+
+  state->f_lower = f_lower;
+  state->f_upper = f_upper;
+
+  if ((f_lower < 0.0 && f_upper < 0.0) || (f_lower > 0.0 && f_upper > 0.0))
+    {
+      GSL_ERROR ("endpoints do not straddle y=0", GSL_EINVAL);
+    }
+
+  return GSL_SUCCESS;
+
+}
+
 static int
 bisection_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper)
 {
@@ -129,6 +149,7 @@ static const gsl_root_fsolver_type bisection_type =
 {"bisection",                           /* name */
  sizeof (bisection_state_t),
  &bisection_init,
+ &bisection_init2,
  &bisection_iterate};
 
 const gsl_root_fsolver_type  * gsl_root_fsolver_bisection = &bisection_type;
diff --git a/roots/brent.c b/roots/brent.c
index 4f0c9cd3..c2a81b09 100644
--- a/roots/brent.c
+++ b/roots/brent.c
@@ -42,6 +42,7 @@ typedef struct
 brent_state_t;
 
 static int brent_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper);
+static int brent_init2 (void * vstate, gsl_function * f, double * root, double x_lower, double f_lower, double x_upper, double f_upper);
 static int brent_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper);
 
 
@@ -69,6 +70,38 @@ brent_init (void * vstate, gsl_function * f, double * root, double x_lower, doub
   state->d = x_upper - x_lower ;
   state->e = x_upper - x_lower ;
 
+  printf("# BRENT INIT with %e %e -- %e %e\n", x_lower, f_lower, x_upper, f_upper);
+
+  if ((f_lower < 0.0 && f_upper < 0.0) || (f_lower > 0.0 && f_upper > 0.0))
+    {
+      GSL_ERROR ("endpoints do not straddle y=0", GSL_EINVAL);
+    }
+
+  return GSL_SUCCESS;
+
+}
+
+static int
+brent_init2 (void * vstate, gsl_function * f, double * root, double x_lower, double f_lower, double x_upper, double f_upper)
+{
+  brent_state_t * state = (brent_state_t *) vstate;
+
+  *root = 0.5 * (x_lower + x_upper) ;
+
+  state->a = x_lower;
+  state->fa = f_lower;
+
+  state->b = x_upper;
+  state->fb = f_upper;
+
+  state->c = x_upper;
+  state->fc = f_upper;
+
+  state->d = x_upper - x_lower ;
+  state->e = x_upper - x_lower ;
+
+  printf("# BRENT INIT with %e %e -- %e %e\n", x_lower, f_lower, x_upper, f_upper);
+
   if ((f_lower < 0.0 && f_upper < 0.0) || (f_lower > 0.0 && f_upper > 0.0))
     {
       GSL_ERROR ("endpoints do not straddle y=0", GSL_EINVAL);
@@ -239,6 +272,7 @@ static const gsl_root_fsolver_type brent_type =
 {"brent",                               /* name */
  sizeof (brent_state_t),
  &brent_init,
+ &brent_init2,
  &brent_iterate};
 
 const gsl_root_fsolver_type  * gsl_root_fsolver_brent = &brent_type;
diff --git a/roots/falsepos.c b/roots/falsepos.c
index f6cda5c8..95080173 100644
--- a/roots/falsepos.c
+++ b/roots/falsepos.c
@@ -52,6 +52,7 @@ typedef struct
 falsepos_state_t;
 
 static int falsepos_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper);
+static int falsepos_init2 (void * vstate, gsl_function * f, double * root, double x_lower, double f_lower, double x_upper, double f_upper);
 static int falsepos_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper);
 
 static int
@@ -78,6 +79,25 @@ falsepos_init (void * vstate, gsl_function * f, double * root, double x_lower, d
 
 }
 
+static int
+falsepos_init2 (void * vstate, gsl_function * f, double * root, double x_lower, double f_lower, double x_upper, double f_upper)
+{
+  falsepos_state_t * state = (falsepos_state_t *) vstate;
+
+  *root = 0.5 * (x_lower + x_upper);
+
+  state->f_lower = f_lower;
+  state->f_upper = f_upper;
+
+  if ((f_lower < 0.0 && f_upper < 0.0) || (f_lower > 0.0 && f_upper > 0.0))
+    {
+      GSL_ERROR ("endpoints do not straddle y=0", GSL_EINVAL);
+    }
+
+  return GSL_SUCCESS;
+
+}
+
 static int
 falsepos_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper)
 {
@@ -173,6 +193,7 @@ static const gsl_root_fsolver_type falsepos_type =
 {"falsepos",                            /* name */
  sizeof (falsepos_state_t),
  &falsepos_init,
+ &falsepos_init2,
  &falsepos_iterate};
 
 const gsl_root_fsolver_type  * gsl_root_fsolver_falsepos = &falsepos_type;
diff --git a/roots/fsolver.c b/roots/fsolver.c
index 8f586bb6..8d4b7ad4 100644
--- a/roots/fsolver.c
+++ b/roots/fsolver.c
@@ -66,6 +66,22 @@ gsl_root_fsolver_set (gsl_root_fsolver * s, gsl_function * f, double x_lower, do
   return (s->type->set) (s->state, s->function, &(s->root), x_lower, x_upper);
 }
 
+int
+gsl_root_fsolver_set2 (gsl_root_fsolver * s, gsl_function * f, double x_lower, double f_lower, double x_upper, double f_upper)
+{
+  if (x_lower > x_upper)
+    {
+      GSL_ERROR ("invalid interval (lower > upper)", GSL_EINVAL);
+    }
+
+  s->function = f;
+  s->root = 0.5 * (x_lower + x_upper);  /* initial estimate */
+  s->x_lower = x_lower;
+  s->x_upper = x_upper;
+
+  return (s->type->set2) (s->state, s->function, &(s->root), x_lower, f_lower, x_upper, f_upper);
+}
+
 int
 gsl_root_fsolver_iterate (gsl_root_fsolver * s)
 {
diff --git a/roots/gsl_roots.h b/roots/gsl_roots.h
index 46e45870..fef4e9e8 100644
--- a/roots/gsl_roots.h
+++ b/roots/gsl_roots.h
@@ -41,6 +41,7 @@ typedef struct
     const char *name;
     size_t size;
     int (*set) (void *state, gsl_function * f, double * root, double x_lower, double x_upper);
+    int (*set2) (void *state, gsl_function * f, double * root, double x_lower, double f_lower, double x_upper, double f_upper);
     int (*iterate) (void *state, gsl_function * f, double * root, double * x_lower, double * x_upper);
   }
 gsl_root_fsolver_type;
@@ -82,6 +83,11 @@ int gsl_root_fsolver_set (gsl_root_fsolver * s,
                           gsl_function * f, 
                           double x_lower, double x_upper);
 
+int gsl_root_fsolver_set2 (gsl_root_fsolver * s,
+                           gsl_function * f,
+                           double x_lower, double f_lower,
+                           double x_upper, double f_upper);
+
 int gsl_root_fsolver_iterate (gsl_root_fsolver * s);
 
 const char * gsl_root_fsolver_name (const gsl_root_fsolver * s);