The initial parameters are n_clusters, n_features, n_cases
 that describe the shape of the problem and n_outer, n_inner
 bounding the loops.

get_int proc (lower_bound in int) ⇒ (result out int)
   get result
   assert result ≥ lower_bound
end proc

n_clusters const = get_int(2)
n_features const = get_int(1)
n_cases    const = get_int(n_clusters)
n_outer    const = get_int(1)
n_inner    const = get_int(1)

 Next we read the feature matrix.

feature var := {{0.0} @@ n_features} @@ n_cases
loop for case_index in [0,n_cases)
    loop for feature_index in [0,n_features)
       get feature[case_index][feature_index]
    end loop
end loop

 Make a random selection of cases to serve as cluster centres.

random_assignment proc (centre out [][]real)
   perm var []int
   loop for i in [0,n_cases)
      perm[i] := i
   end loop
   loop for c in [0,n_clusters)
      n const = n_cases - c                 number of elements to choose from
      i const = n - 1                       last of those elements
      j const = trunc(float(n) * random())  j in [0,n)
      (perm[j], perm[i]) := (perm[i], perm[j])
      centre[c] := feature[perm[i]]
   end loop
end proc    

 The squared Euclidean distance between two feature vectors.

distance_squared proc (x in []real, y in []real) ⇒ (d out real)
   d := 0.0
   loop for i in [0,size x)
      d +:= (x[i] - y[i])**2
   end loop
end proc

repartition proc (
   old_centre in  [][]real,
   new_centre out [][]real,
   group      out []int       case number -> group number
) ⇒ (
   sum_of_squared_distances out real
)
   sum_of_squared_distances := 0.0
   tally var := {0.0} @@ n_clusters
   group := {}
   new_centre := {{0.0} @@ n_features} @@ n_clusters
   loop for case_index in [0,n_cases)
      this_case const = feature[case_index]
      best_distance_squared var := 0.0
      best_cluster          var := -1
      loop for c in [0,n_clusters)
         d const = distance_squared(this_case, old_centre[c])
         if best_cluster < 0 or d < best_distance_squared then
            (best_cluster, best_distance_squared) := (c, d)
         end if
      end if
      assert best_cluster ≥ 0
      group[case_index,) := {best_cluster}
      tally[best_cluster] +:= 1.0
      loop for feature_index in [0,n_features)
         new_centre[best_group, feature_index] +:= this_case[feature_index]
      end loop
      loop for cluster in [0,n_clusters)
         if tally[cluster] > 0.0 then
            loop for feature_index in [0,n_features)
               new_centre[cluster][feature_index] /= tally[cluster]
            end loop
         end if
      end loop
      sum_of_squared_distances +:= best_distance_squared
   end loop
end proc   

outer_best_sum_of_squared_distances var real
outer_best_group var []int

loop for outer in [0,n_outer)
   centre_one var := {{0.0} @@ n_features} @@ n_clusters
   centre_two var := {{0.0} @@ n_features} @@ n_clusters
   d1         var := 0.0
   d2         var := 0.0
   group      var []int

   random_assignment(centre_one)
   loop for inner in [0,n_inner)
      d1 := repartition(centre_one, centre_two, group)
      d2 := repartition(centre_two, centre_one, group)
      while d2 < d1
   end loop
   if outer = 0 or d2 < outer_best_sum_of_squared_distances then
      (outer_best_group, outer_best_sum_of_squared_distances) := (group, d2)
   end if
end loop

print outer_best_group