Implement rudimentary `PowerSemigroup` attribute/operation

[james-d-mitchell]

Something along the lines of:

PowerSemigroup := function(S)
  local SetProduct, PS, table, row, X, Y;

  SetProduct := function(X, Y)
    local result, x, y;

    result := [];
    for x in X do
      for y in Y do
        AddSet(result, x * y);
      od;
    od;
    return result;
  end;
  
  PS := Set(Combinations(AsList(S)));

  table := [];

  for X in PS do
    row := [];
    for Y in PS do
      Add(row, Position(PS, SetProduct(X, Y)));
    od;
    Add(table, row);
  od;
  return SemigroupByMultiplicationTable(table);
end;