function [ x, q0, q1, nl, nf, fx, qf1, qa, qb, qc, qd0, qd1 ] = quad_extrap ( ...
  n, f, x, t, h, v, q0, q1, nl, nf, dmin, ldt, fx, qf1, qa, qb, qc, qd0, qd1 )

%*****************************************************************************80
%
%% quad_extrap() seeks to minimize the scalar function F along a particular curve.
%
%  Discussion:
%
%    The minimizer to be sought is required to lie on a curve defined
%    by Q0, Q1 and X.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    31 July 2016
%
%  Author:
%
%    Original FORTRAN77 version by Richard Brent.
%    This version by John Burkardt.
%
%  Reference:
%
%    Richard Brent,
%    Algorithms for Minimization with Derivatives,
%    Prentice Hall, 1973,
%    Reprinted by Dover, 2002.
%
%  Input:
%
%    integer N, the number of variables.
%
%    real F ( x, n ), the function.
%
%    real X(N), ?
%
%    real T, ?
%
%    real H, ?
%
%    real V(N,N), the matrix of search directions.
%
%    real Q0(N), an auxiliary point used to define
%    a curve through X.
%
%    real Q1(N), an auxiliary point used to define
%    a curve through X.
%
%    integer NL, the number of linear searches.
%
%    integer NF, the number of function evaluations.
%
%    real DMIN, an estimate for the smallest eigenvalue.
%
%    real LDT, the length of the step.
%
%    real FX, the value of F(X,N).
%
%    real QF1, QA, QB, QC, QD0, QD1, ?
%
%  Output:
%
%    real X(N), ?
%
%    real Q0(N), an auxiliary point used to define a curve through X.
%
%    real Q1(N), an auxiliary point used to define a curve through X.
%
%    integer NL, the number of linear searches.
%
%    integer NF, the number of function evaluations.
%
%    real FX, the value of F(X,N).
%
%    real QF1, QA, QB, QC, QD0, QD1, ?
%
  temp = fx;
  fx   = qf1;
  qf1  = temp;

  temp = x;
  x = q1;
  q1 = temp;

  qd1 = norm ( x - q1 );

  l = qd1;
  s = 0.0;

  if ( qd0 <= 0.0 || qd1 <= 0.0 || nl < 3 * n * n )

    fx = qf1;
    qa = 0.0;
    qb = 0.0;
    qc = 1.0;

  else

    jsearch = -1;
    nits = 2;
    value = qf1;
    fk = true;

    [ s, l, value, x, nl, nf, fx, qa, qb, qc ] = minny ( n, jsearch, nits, ...
      s, l, value, fk, f, x, t, h, v, q0, q1, nl, nf, dmin, ldt, ...
      fx, qa, qb, qc, qd0, qd1 );

    qa =                 l * ( l - qd1 )       / ( qd0 + qd1 ) / qd0;
    qb = - ( l + qd0 )     * ( l - qd1 ) / qd1                 / qd0;
    qc =   ( l + qd0 ) * l               / qd1 / ( qd0 + qd1 );

  end

  qd0 = qd1;

  xnew(1:n) = qa * q0(1:n) + qb * x(1:n) + qc * q1(1:n);

  q0(1:n) = x(1:n);
  x(1:n) = xnew(1:n);

  return
end

