``` 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
```