// compute minimum spanning tree via Prim's algorithm
prim := function(A)
  n := Nrows(A);
  V := { 1 }; // start with a tree with a single vertex
  E := []; // empty list of edges
  while #V lt n do
    min := 2^32 + 1;
    // find the shortest edge with one vertex in the tree constructed so far
    // and the other outside of it
    for x in V do
      for y in {1..n} diff V do
        if A[x,y] lt min then
          min := A[x,y];
          e := [x,y];
        end if;
      end for;
    end for;
    // add the new vertex to the tree and the edge to the list of edges
    V := V join {e[2]};
    E := E cat [e];
  end while;
  print "weights of minimal spanning tree:",
        Sort([ A[E[i,1], E[i,2]] : i in [1..#E] ]),
        "total weight", &+[ A[E[i,1], E[i,2]] : i in [1..#E] ];
  return E;
end function;



// compute shortest paths from a to any other point
// if only the shortest path to a single point b is required,
// give b as an additional argument and adjust as indicated
dijkstra := function(A, a) // for paths to all points
// dijkstra := function(A, a, b) // for paths to a single point
  n := Nrows(A);
  S := {}; // set for points which have been dealt with
  d := [ 2^32 : i in [1..n] ]; // array for distances
  p := [ [] : i in [1..n] ]; // array for paths
  d[a] := 0;
  p[a] := [0]; // 0 indicates begin of path
  while #S lt n do // for paths to all points
//  while not b in S do // for paths to a single point
    dmin := 2^32 + 1;
    // find point outside S with minimal distance to a
    for w in {1..n} diff S do
      if d[w] lt dmin then
        v := w;
        dmin := d[v];
      end if;
    end for;
    S := S join {v};
    // check whether going via v to w is shorter than the current best path to w
    for w in {1..n} diff S do
      if d[v] + A[v][w] lt d[w] then
        d[w] := d[v] + A[v][w];
        p[w] := p[v] cat [v];
      end if;
    end for;
  end while;
  return d, p; // for paths to all points
//  return d[b], p[b]; // for paths to a single point
end function;



// compute short cycle via a greedy algorithm, starting with a random vertex
greedy := function(A)
  n := Nrows(A);
  cycle := [ Random([1..n]) ]; // start with a random vertex
  while #cycle lt n do
    // find the closest neighbour not yet in the cycle
    min := 2^31; // arbitrary large number
    for j in [1..n] do
      if not j in cycle and A[cycle[#cycle], j] lt min then
        min := A[cycle[#cycle], j];
        closest := j;
      end if;
    end for;
    cycle := cycle cat [closest]; // append closest neighbour to the cycle
  end while;
  // compute and print the total length of the cycle
  dist := &+[ A[cycle[i], cycle[i+1]] : i in [1..n-1] ] + A[cycle[n], cycle[1]];
  print "greedy algorithm found path of length", dist;
  return cycle;
end function;

