#include <cmath>
#include <cassert>

std::pair<double, double> golden_ratio_search( double f( double x ),
                            double a, double b,
                            double eps_step ) {
    assert( b > a );
    // Calculate these explicitly to maximize precision
    //  - 'phi' is the golden ratio
    //  - 'inv_phi1' is 1/phi
    //  - 'inv_phi2' is 1/phi^2
    double const inv_phi1{ (std::sqrt(5.0) - 1.0)/2.0 };
    double const inv_phi2{ 1 - inv_phi1 };

    double  fa{f( a )};
    double  c1{a + inv_phi2*(b - a)};
    double fc1{f( c1 )};
    double  c2{a + inv_phi1*(b - a)};
    double fc2{f( c2 )};
    double  fb{f( b )};

    while ( (b - a) > eps_step ) {
        if ( fc1 > fc2 ) {
              b = c2;
             fb = fc2;
             c2 = c1;
            fc2 = fc1;
             c1 = a + inv_phi2*(b - a);
            fc1 = f( c1 );
        } else if ( fc1 < fc2 ) {
              a = c1;
             fa = fc1;
             c1 = c2;
            fc1 = fc2;
             c2 = a + inv_phi1*(b - a);
            fc2 = f( c2 );
        } else {
            assert( false );
        }
    }

    // Return the maximum of f(a), f(c ), f(c ) and f(b)
    //                                1      2

    if ( (fa > fc1) && (fa > fc2) && (fa > fb) ) {
        return std::make_pair( a, fa );
    } else if ( (fc1 > fc2) && (fc1 > fb) ) {
        return std::make_pair( c1, fc1 );
    } else if (  fc2 > fb ) {
        return std::make_pair( c2, fc2 );
    } else {
        return std::make_pair( b, fb );
    }
}