function tsp_brute ( filename )

%*****************************************************************************80
%
%% tsp_brute(): brute force solution of small traveling salesman problems.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    30 September 2022
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string FILENAME, the name of the file containing the distance
%    matrix information.  This is an optional input.  If it is not provided,
%    the program will prompt for it.
%
  tic ( )
%
%  Get the distance table.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The program needs to load a file containing\n' );
    fprintf ( 1, '  a city-to-city distance matrix D.\n' );  
    filename = input ( '  Enter the file name, in quotes:  ' );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Reading distance matrix from "%s"\n', filename );
  end

  distance = load ( filename );
%
%  Approve the distance table.
%
  [ m, n ] = size ( distance );

  if ( m ~= n )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'tsp_brute(): Fatal error!\n' );
    fprintf ( 1, '  The distance matrix must be square.\n' );
    fprintf ( 1, '  Your matrix has M = %d, N = %d\n', m, n );
    error ( 'tsp_brute(): Fatal error!\n' );
  end

  v = diag ( distance );
  test = norm ( v );

  if ( 0.0 < test )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'tsp_brute(): Fatal error!\n' );
    fprintf ( 1, '  The distance matrix must have zero diagonal.\n' );
    fprintf ( 1, '  Your matrix has |diag(D)| = %g\n', test );
    error ( 'tsp_brute(): Fatal error!\n' );
  end

  test = norm ( distance - distance' );

  if ( 0.0 < test )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'tsp_brute(): Fatal error!\n' );
    fprintf ( 1, '  The distance matrix must be symmetric.\n' );
    fprintf ( 1, '  Your matrix has |D-D''| = %g\n', test );
    error ( 'tsp_brute(): Fatal error!\n' );
  end
%
%  Print the distance matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The city-to-city distance matrix:\n' )
  fprintf ( 1, '\n' );

  disp ( distance )
%
%  Examine every permutation.
%
  total_max = - inf;
  total_min = inf;
  total_ave = 0.0;

  p = [];
  more = false;
  rank = 0;

  p_min = [];

  paths = 0;

  while ( true )

    [ p, more, rank ] = perm1_next3 ( n, p, more, rank );

    if ( ~ more )
      break
    end

    paths = paths + 1;

    total = path_cost ( n, distance, p );

    total_ave = total_ave + total;

    if ( total_max < total )
      total_max = total;
    end

    if ( total < total_min )
      total_min = total;
      p_min(1:n) = p(1:n);
    end

  end

  total_ave = total_ave / paths;
%
%  Report.
%
  fprintf ( 1, '  A minimal length itinerary:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Step  From  To        Distance\n' );
  fprintf ( 1, '\n' );
  for i1 = 1 : n
    i2 = mod ( i1, n ) + 1;
    fprintf ( 1, '  %4d    %2d  %2d  %14.6g\n', ...
      i2, p_min(i1), p_min(i2), distance ( p_min(i1), p_min(i2) ) );
  end
  fprintf ( 1, '  ----    --  --  --------------\n' );
  fprintf ( 1, '  Total:          %14.6g\n', total_min );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of paths checked = %d\n', paths );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Minimum length = %g\n', total_min );
  fprintf ( 1, '  Average length = %g\n', total_ave );
  fprintf ( 1, '  Maximum length = %g\n', total_max );

  toc ( )

  return
end
function cost = path_cost ( n, distance, p );

%*****************************************************************************80
%
%% path_cost() evaluates the cost of a round trip.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    25 October 2018
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N, the number of cities.
%
%    real DISTANCE(N,N), the city to city distance table.
%
%    integer P(N), a permutation of 1:N, the route.
%
%  Output:
%
%    real COST, the cost of the route.
%
  cost = 0.0;
  i1 = n;
  for i2 = 1 : n
    cost = cost + distance ( p(i1), p(i2) );
    i1 = i2;
  end

  return
end
function [ p, more, rank ] = perm1_next3 ( n, p, more, rank )

%*****************************************************************************80
%
%% perm1_next3() computes permutations of (1,...,N).
%
%  Discussion:
%
%    The routine is initialized by calling with MORE = TRUE, in which case
%    it returns the identity permutation.
%
%    If the routine is called with MORE = FALSE, then the successor of the
%    input permutation is computed.
%
%    Trotter's algorithm is used.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    09 November 2018
%
%  Author:
%
%    Original FORTRAN77 version by Hale Trotter,
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Hale Trotter,
%    PERM, Algorithm 115,
%    Communications of the Association for Computing Machinery,
%    Volume 5, 1962, pages 434-435.
%
%  Input:
%
%    integer N, the number of objects being permuted.
%
%    integer P(N).  This is the previous permutation, 
%    if MORE is TRUE.  But on a startup call, with MORE FALSE, the input value 
%    of P is not needed.
%
%    logical MORE.  On input, should be set to FALSE on an
%    initialization call, and on subsequent inputs should have its output
%    value from the previous call.
%
%    integer RANK, the rank of the current permutation.
%
%  Output:
%
%    integer P(N).  On output, the next permutation; if MORE is FALSE, 
%    then P is the first permutation in the sequence.
%
%    logical MORE.  TRUE if there was a next
%    permutation to produce, or FALSE if there are no more permutations 
%    to produce.
%
%    integer  RANK, the rank of the current permutation.
%
  if ( ~ more )

    p = 1 : n;
    more = true;
    rank = 1;

  else

    n2 = n;
    m2 = rank;
    s = n;

    while ( true )

      q = mod ( m2, n2 );
      t = mod ( m2, 2 * n2 );

      if ( q ~= 0 )
        break
      end

      if ( t == 0 )
        s = s - 1;
      end

      m2 = floor ( m2 / n2 );
      n2 = n2 - 1;

      if ( n2 == 0 )
        p = 1 : n;
        more = false;
        rank = 1;
        break
      end

    end

    if ( n2 ~= 0 )

      if ( q == t )
        s = s - q;
      else
        s = s + q - n2;
      end

      t      = p(s);
      p(s)   = p(s+1);
      p(s+1) = t;

      rank = rank + 1;

    end

  end

  return
end

