function xmin = cubic_interp(x1, x2, f1, f2, g1, g2) 
% call: xmin = cubic_interp(x1, x2, f1, f2, g1, g2) 
% find minimizer of the cubic polynomial Hermite-interpolating a
% function of one variable, at the two points x1 and x2, using the
% function (f1 and f2) and derivative (g1 and g2) data.
% could easily return inf or NaN, so this must be checked after call

% the formula is given by Nocedal and Wright p. 57, but they have
% an error.  Their formula is correct only if x2 > x1, as they
% omitted the sign factor in the definition of gamma.
% Found the correct formula in a 1974 NPL report by Gill and Murray.
% author: M. Overton (overton@cs.nyu.edu), last revised March 2003

eta = g1 + g2 - 3*(f1 - f2)/(x1 - x2);    % called d1 in Nocedal and Wright
gamma = sign(x2-x1)*sqrt(eta^2 - g1*g2);  % called d2 in Nocedal and Wright
xmin = x2 - (x2 - x1)*(g2 + gamma - eta)/(g2 - g1 + 2*gamma);
   


