function [ x_new, y_new ] = diophantine_solution_minimize ( a, b, x, y )

%*****************************************************************************80
%
%% DIOPHANTINE_SOLUTION_MINIMIZE seeks a minimal solution of a Diophantine equation.
%
%  Discussion:
%
%    Given a solution (X,Y) of a Diophantine equation:
%
%      A * X + B * Y = C.
%
%    then there are an infinite family of solutions of the form
%
%      ( X(i), Y(i) ) = ( X + i * B, Y - i * A )
%
%    An integral solution of minimal Euclidean norm can be found by
%    tentatively moving along the vectors (B,-A) and (-B,A) one step
%    at a time.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Eric Weisstein, editor,
%    CRC Concise Encylopedia of Mathematics,
%    CRC Press, 1998, page 446.
%
%  Parameters:
%
%    Input, integer A, B, the coefficients of the Diophantine equation.
%    A and B are assumed to be relatively prime.
%
%    Input, integer X, Y, on input, a solution of the Diophantine equation.  
%
%    Output, integer X_NEW, Y_NEW, a solution of minimal Euclidean norm.
%
  x_new = x;
  y_new = y;
%
%  Compute the minimum for T real, and then look nearby.
%
  t = ( - b * x_new + a * y_new ) / ( a * a + b * b );

  x_new = x_new + round ( t ) * b;
  y_new = y_new - round ( t ) * a;
%
%  Look nearby.
%
  norm = x_new * x_new + y_new * y_new;

  while ( 1 )

    x2 = x_new + b;
    y2 = y_new - a;

    norm2 = x2 * x2 + y2 * y2;

    if ( norm <= norm2 )
      break
    end

    x_new = x2;
    y_new = y2;
    norm = norm2;

  end

  while ( 1 )

    x2 = x_new - b;
    y2 = y_new + a;

    norm2 = x2 * x2 + y2 * y2;

    if ( norm <= norm2 )
      break
    end

    x_new = x2;
    y_new = y2;
    norm = norm2;

  end

  return
end