function EC = jequiv_classes(S)
% for all a,b in S, a J b   iff SaS = SbS
% brute force

[n,n] = size(S);
EQ = zeros(n);
for i=1:n
  for j=i+1:n
    Si = Sab(S,i);
    Sj = Sab(S,j);
    EQ(i,j) = seteq(Si,Sj);
  end
end

EQ = EQ+EQ'+eye(n);
X = 1:n;
EC={};
while ~isempty(X)
    x = X(1);
    Ex = find(EQ(x,:));
    EC{end+1} = Ex;
    X = setdiff(X,Ex);
end
return



function Si = Sab(S,a)
[n,n] = size(S);
Si = [];
for i=1:n
  for j=1:n
    Si(end+1) = S(S(i,a),j);
  end
end
return