# include <stdlib.h>
# include <stdio.h>

# include "bisection_rc.h"

/******************************************************************************/

double bisection_rc ( double *a, double *b, double fx, int *job )

/******************************************************************************/
/*
  Purpose:

    bisection_rc() seeks a zero of f(x) in a change of sign interval.

  Discussion:

    The bisection method is used.

    This routine uses reverse communication, so that the function is always
    evaluated in the calling program.

    On the first call, the user sets JOB = 0, and the values of A and B.
    Thereafter, the user checks the returned value of JOB and follows 
    directions.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    14 January 2013

  Author:

    John Burkardt

  Parameters:

    Input/output, double *A, *B, the endpoints of the change of 
    sign interval.  These values will change from call to call as the
    interval size is reduced.

    Input, double FX, the function value at the point X returned on
    the previous call.  On the first call, FX is ignored.

    Input/output, int *JOB, a communication flag.
    The user sets JOB to 0 before the first call.  Thereafter, the program
    controls setting the value of JOB, whose output values mean:

    Output, double BISECTION_RC, a point X at which the function is to 
    be evaluated.
*/
{
  static double fa;
  static double fb;
  static int state;
  double x;
  static double x_local;

  if ( *job == 0 )
  {
    fa = 0.0;
    fb = 0.0;
    state = 1;
    x = *a;
    *job = 1;
  }
  else if ( state == 1 )
  {
    if ( fx == 0.0 )
    {
      b = a;
      state = 3;
      return *b;
    }
    fa = fx;
    x = *b;
    state = 2;
  }
  else if ( state == 2 )
  {
    if ( fx == 0.0 )
    {
      a = b;
      state = 3;
      return *b;
    }
    fb = fx;

    if ( r8_sign ( fa ) == r8_sign ( fb ) )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "bisection_rc(): Fatal error!\n" );
      fprintf ( stderr, "  F(A) and F(B) have the same sign.\n" );
      exit ( 1 );
    }

    x = ( *a + *b ) / 2.0;
    state = 3;
  }
  else
  {
    if ( r8_sign ( fx ) == r8_sign ( fa ) )
    {
      *a = x_local;
      fa = fx;
    }
    else
    {
      *b = x_local;
      fb = fx;
    }
    x = ( *a + *b ) / 2.0;
    state = 3;
  }

  x_local = x;

  return x;
}
/******************************************************************************/

double r8_sign ( double x )

/******************************************************************************/
/*
  Purpose:

    R8_SIGN returns the sign of an R8.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    08 May 2006

  Author:

    John Burkardt

  Parameters:

    Input, double X, the number whose sign is desired.

    Output, double R8_SIGN, the sign of X.
*/
{
  double value;

  if ( x < 0.0 )
  {
    value = - 1.0;
  }
  else
  {
    value = + 1.0;
  }
  return value;
}