# include <cmath>
# include <cstdlib>

using namespace std;

# include "golden_section.hpp"

//****************************************************************************80

void golden_section ( double f ( double x ), int n, double x_tol, double &a, 
  double &b, int &it )

//****************************************************************************80
//
//  Purpose:
//
//    golden_section() seeks a minimizer of a unimodal function in [a,b].
//
//  Discussion:
//
//    A unimodal function f(x) in [a,b] has a minimizer a <= c <= b such that
//    f(x) strictly decreases between a and c, and strictly increases between
//    c and b.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 February 2026
//
//  Author:
//
//    Original MATLAB code by J Nathan Kutz.
//    This version by John Burkardt.
//
//  Reference:
//
//    J Nathan Kutz,
//    Data Driven Modeling and Scientific Computation,
//    Oxford University Press, 2013,
//    ISBN: 978-0-19-966034-6.
//
//  Input:
//
//    double f ( double x ): the function to be minimized.
//
//    int n: the maximum number of iterations allowed.
//    n must be at least 1.
//
//    double x_tol: a tolerance for the width of the x interval.
//
//    double &a, &b: the left and right endpoints.
//
//  Output:
//
//    double &a, &b: the current interval containing the minimizer.
//
//    int &it: the number of steps taken.
//
{
  double f1;
  double f2;
  double g;
  double x1;
  double x2;

  for ( it = 0; it <= n; it++ )
  {
    if ( it == 0 )
    {
      g = ( - 1.0 + sqrt ( 5.0 ) ) / 2.0;
      x1 =         g   * a + ( 1.0 - g ) * b;
      x2 = ( 1.0 - g ) * a +         g   * b;
      f1 = f ( x1 );
      f2 = f ( x2 );
    }
    else
    {
      if ( f1 < f2 )
      {
        b = x2;
        x2 = x1;
        f2 = f1;
        x1 = g * a + ( 1.0 - g ) * b;
        f1 = f ( x1 );
      }
      else
      {
        a = x1;
        x1 = x2;
        f1 = f2;
        x2 = ( 1.0 - g ) * a + g * b;
        f2 = f ( x2 );
      }
    }

    if ( fabs ( b - a ) <= x_tol )
    {
      break;
    }
  }

  return;
}