# include <float.h>
# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <time.h>

# include "root_rc.h"

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

double root_rc ( double x, double fx, double *ferr, double *xerr, double q[9] )

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

    ROOT_RC solves a single nonlinear equation using reverse communication.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    22 January 2013

  Author:

    Original FORTRAN77 version by Gaston Gonnet.
    C version by John Burkardt.

  Reference:
     
    Gaston Gonnet,
    On the Structure of Zero Finders,
    BIT Numerical Mathematics,
    Volume 17, Number 2, June 1977, pages 170-183.

  Parameters:

    Input, double X, an estimate for the root.  On the first
    call, this must be a value chosen by the user.  Thereafter, it may
    be a value chosen by the user, or the value of ROOT returned on the
    previous call to the function.

    Input, double FX, the value of the function at X.

    Output, double *FERR, the smallest value of F encountered.

    Output, double *XERR, the width of the change-in-sign interval,
    if one was encountered.

    Input/output, double Q[9], storage needed by the function.
    Before the first call, the user must set Q(1) to 0.

    Output, double ROOT_RC, an improved estimate for the root.
*/
{
  double d;
  double decr;
  int i;
  double p;
  double r;
  double u;
  double v;
  double w;
  double xnew;
  double z;
/*
  If we found an exact zero, there is nothing more to do.
*/
  if ( fx == 0.0 )
  {
    *ferr = 0.0;
    *xerr = 0.0;
    xnew = x;
    return xnew;
  }

  *ferr = fabs ( fx );
/*
  If this is the first time, initialize, estimate the first root, and exit.
*/
  if ( q[0] == 0.0 )
  {
    q[0] = fx;
    q[1] = x;
    for ( i = 2; i < 9; i++ )
    {
      q[i] = 0.0;
    }
    xnew = x + fx;
    *xerr = HUGE_VAL;
    return xnew;
  }
/*
  This is not the first call.
*/
  q[8] = q[8] + 1.0;
/*
  Check for too many iterations.
*/
  if ( 80.0 < q[8] )
  {
    printf ( "\n" );
    printf ( "ROOT_RC - Fatal error!\n" );
    printf ( "  Number of iterations = %d\n", ( int ) q[8] );
    exit ( 1 );
  }
/*
  Check for a repeated X value.
*/
  if ( ( 2.0 <= q[8] && x == q[3] ) || x == q[1] )
  {
    printf ( "\n" );
    printf ( "ROOT_RC - Fatal error!\n" );
    printf ( "  Value of X has been input before.\n" );
    exit ( 1 );
  }
/*
  Push X -> A -> B -> C
*/
  for ( i = 5; 2 <= i; i-- )
  {
    q[i] = q[i-2];
  }
  q[0] = fx;
  q[1] = x;
/*
  If we have a change-in-sign interval, store the opposite value.
*/
  if ( r8_sign ( q[0] ) != r8_sign ( q[2] ) )
  {
    q[6] = q[2];
    q[7] = q[3];
  }
/*
  Calculate XERR.
*/
  if ( q[6] != 0.0 )
  {
    *xerr = fabs ( q[7] - q[1] );
  }
  else
  {
    *xerr = HUGE_VAL;
  }
/*
  If more than 30 iterations, and we have change-in-sign interval, bisect.
*/
  if ( 30.0 < q[8] && q[6] != 0.0 )
  {
    xnew = q[1] + ( q[7] - q[1] ) / 2.0;
    return xnew;
  }

  v = ( q[2] - q[0] ) / ( q[3] - q[1] );
/*
  If 3 or more points, try Muller.
*/
  if ( q[4] != 0.0 )
  {
    u = ( q[4] - q[2] ) / ( q[5] - q[3] );
    w = q[3] - q[1];
    z = ( q[5] - q[1] ) / w;
    r = ( z + 1.0 ) * v - u;

    if ( r != 0.0 )
    {
      p = 2.0 * z * q[0] / r;
      d = 2.0 * p / ( w * r ) * ( v - u );
      if ( -1.0 <= d )
      {
        xnew = q[1] - p / ( 1.0 + sqrt ( 1.0 + d ) );
        if ( q[6] == 0.0 || 
             ( q[1] < xnew && xnew < q[7] ) || 
             ( q[7] < xnew && xnew < q[1] ) )
        {
          return xnew;
        }
      }
    }
  }
/*
  Try the secant step.
*/
  if ( q[0] != q[2] || q[6] == 0.0 )
  {
    if ( q[0] == q[2] )
    {
      printf ( "\n" );
      printf ( "ROOT_RC - Fatal error!\n" );
      printf ( "  Cannot apply any method.\n" );
      exit ( 1 );
    }
    decr = q[0] / v;
    if ( fabs ( decr ) * 4.6E+18 < fabs ( q[1] ) )
    {
      decr = 1.74E-18 * fabs ( q[1] ) * r8_sign ( decr );
    }
    xnew = q[1] - decr;
    if ( q[6] == 0.0 || 
        ( q[1] < xnew && xnew < q[7] ) || 
        ( q[7] < xnew && xnew < q[1] ) )
    {
      return xnew;
    }
  }
/*
  Apply bisection.
*/
  xnew = q[1] + ( q[7] - q[1] ) / 2.0;

  return xnew;
}
/******************************************************************************/

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;
}
/******************************************************************************/

void timestamp ( )

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

    TIMESTAMP prints the current YMDHMS date as a time stamp.

  Example:

    31 May 2001 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    24 September 2003

  Author:

    John Burkardt

  Parameters:

    None
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  fprintf ( stdout, "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}