/* Bracketing function's root. The function is supposed to have unlimited 
   domain and be continuous.

   int BracketRoot(double x0,double *a,double *b,double d0,
                    double di, double dmax,double (*fun)(double));
   Parameters:
    x0 - initial guess on input;
    a  - left bound on output;
    b  - right bound on output;
    d0 - initial interval of hunting;
    di - interval increment (geometric progression multiplier);
    dmax - maximal increment;
    fun - pointer to the function.
   Returns: 
    1 - if a root is bracketed;
    0 - on failure
*/

int BracketRoot(double x0,double *a,double *b,double d0,
                 double di, double dmax,double (*fun)(double)) {
  double fa,fb,f0;
  /* get initial function guess, initial a,b,fa,fb */
  f0=(*fun)(x0); *a=x0-d0; *b=x0+d0; fa=(*fun)(*a); fb=(*fun)(*b);
  /* while the increased search interval is less than maximal, 
     process cycle */
  while((d0*=di)<dmax) {
    /* check up the bracketing success. Case f0>0. */
    if(f0>=0.) {
      if(fa<0.) {*b=x0;return(1);}
      if(fb<0.) {*a=x0;return(1);}
      /* else, compare fa and fb, choose the direction of search. The 
         right search case. */
      if(fa>fb) {*a=x0; x0=(*b); *b+=d0; fa=f0; f0=fb; fb=(*fun)(*b);}
      /* the left search case */
      else if(fa<fb) {*b=x0; x0=(*a); *a-=d0; fb=f0; f0=fa; fa=(*fun)(*a);}
      /* both sides search */
      else {*a-=d0; *b+=d0; fa=(*fun(*a);fb=(*fun)(*b);}
    }
    /* Analogically, case when f0<0 */
    else if(f0<0.) {
      if(fa>=0.) {*b=x0;return(1);}
      else if(fb>=0.) {*a=x0;return(1);}
      /* else, compare fa and fb, choose the direction of search. The 
         right search case. */
      if(fa<fb) {*a=x0; x0=(*b); *b+=d0; fa=f0; f0=fb; fb=(*fun)(*b);}
      /* the left search case */
      else if(fa>fb) {*b=x0; x0=(*a); *a-=d0; fb=f0; f0=fa; fa=(*fun)(*a);}
      /* both sides search */
      else {*a-=d0; *b+=d0; fa=(*fun(*a);fb=(*fun)(*b);}
    }
  }
  /* if we get there, the search failed */
  return(0);
}