subroutine tsp_anneal ( n, x, cost, order )

!*****************************************************************************80
!
!! tsp_anneal() solves a traveling salesman problem by simulated annealing.
!
!  Discussion:
!
!    The code seeks the shortest round-trip to a number of cities whose 
!    coordinates are known.  This is the traveling salesman problem (TSP).
!    The simulated annealing method is used, which is not guaranteed to find
!    the true best path.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    19 June 2026
!
!  Author:
!
!    Original Fortran77 version by Press, Flannery, Teukolsky, Vetterling.
!    This version by John Burkardt.
!
!  Reference:
!
!    William Press, Brian Flannery, Saul Teukolsky, William Vetterling,
!    Numerical Recipes in FORTRAN: The Art of Scientific Computing,
!    Third Edition,
!    Cambridge University Press, 2007,
!    ISBN: 978-0-521-88068-8,
!    LC: QA297.N866.
!
!  Input:
!
!    integer n: the number of cities.
!
!    real ( kind = rk8 ) x(n,2): the coordinates of the cities.
!
!  Output:
!
!    real ( kind = rk8 ) cost: the cost of the best path found.
!
!    integer order(n): the best path found.
! 
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  logical accept
  real ( kind = rk8 ) cost
  real ( kind = rk8 ) de
  integer j
  integer k
  integer nc(6)
  integer nlimit
  integer nover
  integer ns
  integer nsucc
  integer order(n)
  real ( kind = rk8 ) r
  real ( kind = rk8 ) t
  real ( kind = rk8 ), parameter :: tfactr = 0.9
  real ( kind = rk8 ) x(n,2)
 
  nover = 100 * n
  nlimit = 10 * n
  t = 0.5
!
!  Initialize the ordering.
!
  call perm1_random ( n, order )
!
!  Compute the cost of the initial path.
!
  call tsp_tour_cost ( n, x, order, cost )

  write ( *, '(a,g14.6)' ) '  Cost of initial random tour is ', cost
!
!  Generate samples at 100 decreasing temperatures.
!
  do j = 1, 100

    nsucc = 0
!
!  At each temperature, make NOVER tries.
!
    do k = 1, nover
!
!  Randomly choose two cities nc(1) and nc(2).
!
      do

        call random_number ( harvest = r )
        nc(1) = 1 + int ( n * r )

        call random_number ( harvest = r )
        nc(2) = 1 + int ( ( n - 1 ) * r )

        if ( nc(1) <= nc(2) ) then
          nc(2) = nc(2) + 1
        end if

        ns = 1 + mod ( nc(1) - nc(2) + n - 1, n )

        if ( 3 <= ns ) then
          exit
        end if

      end do
!
!  Randomly consider either a transposition or a reversal.
!
      call random_number ( harvest = r )

      if ( r < 0.5 ) then

        call random_number ( harvest = r )
        nc(3) = nc(2) + int ( abs ( ns - 2 ) * r ) + 1
        nc(3) = 1 + mod ( nc(3) - 1, n )

        call transpose_cost ( x, order, n, nc, de )

        call cost_acceptance ( de, t, accept )

        if ( accept ) then
          nsucc = nsucc + 1
          call transpose_do ( order, n, nc )
          call tsp_tour_cost ( n, x, order, cost )
        end if

      else

        call reversal_cost ( x, order, n, nc, de )

        call cost_acceptance ( de, t, accept )

        if ( accept ) then
          nsucc = nsucc + 1
          call reversal_do ( order, n, nc )
          call tsp_tour_cost ( n, x, order, cost )
        end if

      end if

      if ( nlimit <= nsucc ) then
        exit
      end if

    end do

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'tsp_anneal():'
    write ( *, '(a,f10.6)' ) '  Temperature = ', t
    write ( *, '(a,f10.6)' ) '  Tour cost  = ', cost
    write ( *, '(a,i6)' ) '  Successful moves: ', nsucc

    if ( nsucc == 0 ) then
      return
    end if
!
!  Lower the temperature for the next trial.
!
    t = t * tfactr

  end do

  return
end
subroutine cost_acceptance ( de, t, accept )

!*****************************************************************************80
!
!! cost_acceptance() chooses whether to accept a proposed alternate solution.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    20 June 2026
!
!  Author:
!
!    Original Fortran77 version by Press, Flannery, Teukolsky, Vetterling.
!    This version by John Burkardt.
!
!  Reference:
!
!    William Press, Brian Flannery, Saul Teukolsky, William Vetterling,
!    Numerical Recipes in FORTRAN: The Art of Scientific Computing,
!    Third Edition,
!    Cambridge University Press, 2007,
!    ISBN: 978-0-521-88068-8,
!    LC: QA297.N866.
!
!  Input:
!
!    real ( kind = rk8 ) de: the cost of a change in the ordering.
!
!    real ( kind = rk8 ) t: the current "temperature".
!
!  Output:
!
!    logical accept: TRUE if the cost change is acceptable.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  logical accept
  real ( kind = rk8 ) de
  real ( kind = rk8 ) r
  real ( kind = rk8 ) t

  if ( de < 0.0E+00 ) then
    accept = .true.
  else
    call random_number ( harvest = r )
    accept = ( r < exp ( - de / t ) )
  end if

  return
end
function i4_uniform_ab ( a, b )

!*****************************************************************************80
!
!! i4_uniform_ab() returns a scaled pseudorandom I4 between A and B.
!
!  Discussion:
!
!    An I4 is an integer value.
!
!    The pseudorandom number will be scaled to be uniformly distributed
!    between A and B.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    02 October 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Paul Bratley, Bennett Fox, Linus Schrage,
!    A Guide to Simulation,
!    Second Edition,
!    Springer, 1987,
!    ISBN: 0387964673,
!    LC: QA76.9.C65.B73.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, December 1986, pages 362-376.
!
!    Pierre L'Ecuyer,
!    Random Number Generation,
!    in Handbook of Simulation,
!    edited by Jerry Banks,
!    Wiley, 1998,
!    ISBN: 0471134031,
!    LC: T57.62.H37.
!
!    Peter Lewis, Allen Goodman, James Miller,
!    A Pseudo-Random Number Generator for the System/360,
!    IBM Systems Journal,
!    Volume 8, Number 2, 1969, pages 136-143.
!
!  Input:
!
!    integer A, B, the limits of the interval.
!
!  Output:
!
!    integer I4_UNIFORM_AB, a number between A and B.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer a
  integer b
  integer, parameter :: i4_huge = 2147483647
  integer i4_uniform_ab
  real r
  integer value

  call random_number ( harvest = r )
!
!  Scale R to lie between A-0.5 and B+0.5.
!
  r = ( 1.0E+00 - r ) * ( real ( min ( a, b ) ) - 0.5E+00 ) & 
    +             r   * ( real ( max ( a, b ) ) + 0.5E+00 )
!
!  Use rounding to convert R to an integer between A and B.
!
  value = nint ( r )

  value = max ( value, min ( a, b ) )
  value = min ( value, max ( a, b ) )

  i4_uniform_ab = value

  return
end
subroutine perm1_random ( n, p )

!*****************************************************************************80
!
!! perm1_random() selects a random permutation of (1,...,N).
!
!  Discussion:
!
!    The algorithm is known as the Fisher-Yates or Knuth shuffle.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    12 May 2002
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Albert Nijenhuis, Herbert Wilf,
!    Combinatorial Algorithms for Computers and Calculators,
!    Second Edition,
!    Academic Press, 1978,
!    ISBN: 0-12-519260-6,
!    LC: QA164.N54.
!
!  Input:
!
!    integer N, the number of objects to be permuted.
!
!  Output:
!
!    integer P(N), a permutation of ( 1, 2, ..., N ),
!    in standard index form.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  integer i
  integer i4_uniform_ab
  integer j
  integer p(n)
  integer t

  do i = 1, n
    p(i) = i
  end do

  do i = 1, n - 1
    j = i4_uniform_ab ( i, n )
    t    = p(i)
    p(i) = p(j)
    p(j) = t
  end do

  return
end
subroutine reversal_cost ( x, order, n, nc, de )

!*****************************************************************************80
!
!! reversal_cost() evaluates the cost change for a reversal.
!
!  Discussion:
!
!    The current tour includes a segment starting with city 1 and going to
!    city 2.  We plan to reverse this segment.  
!
!    Thus, initially, we have something like:
!
!    --3--1--a--...--z--2--4--
!
!    and after reversal:
!
!    --3--2--z--...--a--1--4--
!
!    The actual identifiers of the cities 1, 2, 3, and 4 are stored
!    in the array nc().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    20 June 2026
!
!  Author:
!
!    Original Fortran77 version by Press, Flannery, Teukolsky, Vetterling.
!    This version by John Burkardt.
!
!  Reference:
!
!    William Press, Brian Flannery, Saul Teukolsky, William Vetterling,
!    Numerical Recipes in FORTRAN: The Art of Scientific Computing,
!    Third Edition,
!    Cambridge University Press, 2007,
!    ISBN: 978-0-521-88068-8,
!    LC: QA297.N866.
!
!  Input:
!
!    real ( kind = rk8 ) x(n,2): the coordinates of the cities.
!
!    integer order(n): the order of visiting the cities.
!
!    integer n: the number of cities.
!
!    integer nc(4): information describing the reversal.
!
!  Output:
!
!    real ( kind = rk8 ) de: the cost of the reversal.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) de
  integer i
  integer j
  integer nc(4)
  integer order(n)
  real ( kind = rk8 ) x(n,2)
  real ( kind = rk8 ) xx(4,2)
!
!  Determine nc(3) and nc(4), the cities initially 1 before and 1 after
!  cities nc(1) and nc(2).
!
  nc(3) = 1 + mod ( ( nc(1) + n - 2 ), n )
  nc(4) = 1 + mod (   nc(2),               n )
!
!  For convenience, make a local copy of the coordinates.
!
  do j = 1, 4
    i = order ( nc(j) )
    xx(j,:) = x(i,:)
  end do
!
!  Compute the total length change by subtracting the deleted lengths, 
!  and adding the new ones.
!
  de = - norm2 ( xx(1,:) - xx(3,:) ) &
       - norm2 ( xx(2,:) - xx(4,:) ) &
       + norm2 ( xx(1,:) - xx(4,:) ) &
       + norm2 ( xx(2,:) - xx(3,:) )

  return
end
subroutine reversal_do ( order, n, nc )

!*****************************************************************************80
!
!! reversal_do() carries out a reversal for the annealing algorithm.
!
!  Discussion:
!
!    The current tour includes a segment starting with city 1 and going to
!    city 2.  We plan to reverse this segment.  
!
!    Thus, initially, we have something like:
!
!    --3--1--a--...--z--2--4--
!
!    and after reversal:
!
!    --3--2--z--...--a--1--4--
!
!    The actual identifiers of the cities 1, 2, 3, and 4 are stored
!    in the array nc.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    20 June 2026
!
!  Author:
!
!    Original Fortran77 version by Press, Flannery, Teukolsky, Vetterling.
!    This version by John Burkardt.
!
!  Reference:
!
!    William Press, Brian Flannery, Saul Teukolsky, William Vetterling,
!    Numerical Recipes in FORTRAN: The Art of Scientific Computing,
!    Third Edition,
!    Cambridge University Press, 2007,
!    ISBN: 978-0-521-88068-8,
!    LC: QA297.N866.
!
!  Input:
!
!    integer order(n): the order of visiting the cities.
!
!    integer n: the number of cities.
!
!    integer nc(4): information describing the reversal.
!
!  Output:
!
!    integer order(n): the updated order of visiting the cities.
!
  implicit none

  integer n

  integer i1
  integer i2
  integer nc(4)
  integer order(n)

  i1 = 1 + mod ( ( nc(1) - 1 ),     n )
  i2 = 1 + mod ( ( nc(2) - 1 + n ), n )

  order(i1:i2) = order(i2:i1:-1)

  return
end
subroutine transpose_cost ( x, order, n, nc, de )

!*****************************************************************************80
!
!! transpose_cost() determines the cost of a transposition for the annealing algorithm.
!
!  Discussion:
!
!    The current tour includes a segment starting with city 1 and going to
!    city 2.  We plan to remove this segment from the current tour, and
!    insert it between the consecutive cities numbered 3 and 4.  Note that
!    city 5 just precedes city 1, and city 6 just follows city 2.  
!
!    Thus, initially, we have something like:
!
!    --5--1--a--...--z--2--6--     --3--4--
!
!    and after transposition:
!
!    --5--6--                      --3--1--a--...--z--2--4--
!
!    The actual identifiers of the cities 1, 2, 3, 4, 5, and 6 are stored
!    in the array nc().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    20 June 2026
!
!  Author:
!
!    Original Fortran77 version by Press, Flannery, Teukolsky, Vetterling.
!    This version by John Burkardt.
!
!  Reference:
!
!    William Press, Brian Flannery, Saul Teukolsky, William Vetterling,
!    Numerical Recipes in FORTRAN: The Art of Scientific Computing,
!    Third Edition,
!    Cambridge University Press, 2007,
!    ISBN: 978-0-521-88068-8,
!    LC: QA297.N866.
!
!  Input:
!
!    real ( kind = rk8 ) x(n,2): the coordinates of the cities.
!
!    integer order(n): the order of visiting the cities.
!
!    integer n: the number of cities.
!
!    integer nc(6): information describing the transposition.
!
!  Output:
!
!    real ( kind = rk8 ) de: the cost of the transposition.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) de
  integer i
  integer j
  integer nc(6)
  integer order(n)
  real ( kind = rk8 ) x(n,2)
  real ( kind = rk8 ) xx(6,2)

  nc(4) = 1 + mod (   nc(3),               n )
  nc(5) = 1 + mod ( ( nc(1) + n - 2 ), n )
  nc(6) = 1 + mod (   nc(2),               n )

  do j = 1, 6
    i = order ( nc(j) )
    xx(j,1:2) = x(i,1:2)
  end do
!
!  Subtract the deleted lengths, and add the new ones.
!
  de = - norm2 ( xx(2,:) - xx(6,:) ) &
       - norm2 ( xx(1,:) - xx(5,:) ) &
       - norm2 ( xx(3,:) - xx(4,:) ) &
       + norm2 ( xx(1,:) - xx(3,:) ) &
       + norm2 ( xx(2,:) - xx(4,:) ) &
       + norm2 ( xx(5,:) - xx(6,:) ) 

  return
end
subroutine transpose_do ( order, n, nc )

!*****************************************************************************80
!
!! transpose_do() carries out a transposition for the annealing algorithm.
!
!  Discussion:
!
!    The current tour includes a segment starting with city 1 and going to
!    city 2.  We plan to remove this segment from the current tour, and
!    insert it between the consecutive cities numbered 3 and 4.  Note that
!    city 5 just precedes city 1, and city 6 just follows city 2.  
!
!    Thus, initially, we have something like:
!
!    --5--1--a--...--z--2--6--     --3--4--
!
!    and after transposition:
!
!    --5--6--                      --3--1--a--...--z--2--4--
!
!    The actual identifiers of the cities 1, 2, 3, 4, 5, and 6 are stored
!    in the array nc().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    20 June 2026
!
!  Author:
!
!    Original Fortran77 version by Press, Flannery, Teukolsky, Vetterling.
!    This version by John Burkardt.
!
!  Reference:
!
!    William Press, Brian Flannery, Saul Teukolsky, William Vetterling,
!    Numerical Recipes in FORTRAN: The Art of Scientific Computing,
!    Third Edition,
!    Cambridge University Press, 2007,
!    ISBN: 978-0-521-88068-8,
!    LC: QA297.N866.
!
!  Input:
!
!    integer order(n): the order of visiting the cities.
!
!    integer n: the number of cities.
!
!    integer nc(6): information describing the transposition.
!
!  Output:
!
!    integer order(n): the revised order of visiting the cities.
!
  implicit none

  integer n

  integer i
  integer j
  integer jj
  integer jorder(n)
  integer m
  integer nc(6)
  integer order(n)

  i = 1
!
!  Cities nc1 to nc2 (=nc4-1).
!
  m = 1 + mod ( ( nc(2) - nc(1) + n ), n )
  do j = 1, m
    jj = 1 + mod ( ( j + nc(1) - 2 ), n )
    jorder(i) = order(jj)
    i = i + 1
  end do
!
!  Cities nc4 to nc5 (=nc6-1)
!
  m = 1 + mod ( ( nc(5) - nc(4) + n ), n )
  do j = 1, m
    jj = 1 + mod ( ( j + nc(4) - 2 ), n )
    jorder(i) = order(jj)
    i = i + 1
  end do
!
!  Cities nc6 to nc3 (=nc1-1)
!
  m = 1 + mod ( ( nc(3) - nc(6) + n ), n )
  do j = 1, m
    jj = 1 + mod ( ( j + nc(6) - 2 ), n )
    jorder(i) = order(jj)
    i = i + 1
  end do
!
!  Replace the old order.
!
  order(1:n) = jorder(1:n)

  return
end
subroutine tsp_tour_cost ( n, x, order, cost )

!*****************************************************************************80
!
!! tsp_tour_cost() evaluates the cost of a round trip.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    19 June 2026
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer n: the number of cities.
!
!    real x(n,2): the city locations.
!
!    integer order(n): a permutation of 1:N, the route.
!
!  Output:
!
!    real cost: the cost of the route.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) cost
  integer from
  integer order(n)
  integer to
  real ( kind = rk8 ) x(n,2)

  cost = 0.0D+00
  from = n
  do to = 1, n
    cost = cost + sqrt ( ( x(order(from),1) - x(order(to),1) )**2 &
                       + ( x(order(from),2) - x(order(to),2) )**2 )
    from = to
  end do

  return
end