Add UniformPartialPerm

Author: james-d-mitchell

PackageInfo.g

@@ -384,7 +384,7 @@ PackageDoc := rec(
 
 Dependencies := rec(
   GAP := ">=4.12.1",
-  NeededOtherPackages := [["datastructures", ">=0.2.5"],
+  NeededOtherPackages := [["datastructures", ">=0.3.0"],
                           ["digraphs", ">=1.6.2"],
                           ["genss", ">=1.6.5"],
                           ["images", ">=1.3.1"],

gap/elements/pperm.gd

@@ -9,3 +9,4 @@
 ##
 
 DeclareAttribute("CyclesOfPartialPerm", IsPartialPerm);
+DeclareOperation("UniformRandomPartialPerm", [IsPosInt]);

gap/elements/pperm.gi

@@ -12,3 +12,70 @@ InstallMethod(IndexPeriodOfSemigroupElement, "for a partial perm",
 
 InstallMethod(CyclesOfPartialPerm, "for a partial perm", [IsPartialPerm],
 f -> DigraphAllSimpleCircuits(AsDigraph(f)));
+
+DeclareGlobalVariable("_RANDOM_PPERM_CACHED_VALUES");
+MakeReadWriteGlobal("_RANDOM_PPERM_CACHED_VALUES");
+_RANDOM_PPERM_CACHED_VALUES := WeakPointerObj([]);
+
+InstallMethod(UniformRandomPartialPerm, "for a positive integer", [IsPosInt],
+function(n)
+  local cache, num_bijections, r, im, unused_vals, num_undefined,
+  num_each_defined, pos, i, j;
+
+  cache := Immutable(_RANDOM_PPERM_CACHED_VALUES);
+
+  if not IsBound(cache[n + 1]) then
+    # TODO check that all the values in cache are bound
+    cache := [];
+    cache[1] := [1];
+    for i in [2 .. n + 1] do
+      cache[i] := [1];
+      for j in [2 .. i] do
+        cache[i][j] := cache[i][j - 1] + (i - 1) * cache[i - 1][j - 1];
+      od;
+    od;
+  fi;
+
+  # The number of partial perms from a set of size i to a set of size j.
+  # num_bijections := function(i, j)
+  #   local tmp;
+  #   if j < i then
+  #     tmp := i;
+  #     i := j;
+  #     j := tmp;
+  #   fi;
+  #   return Sum([0 .. i], r -> Binomial(i, r) * Binomial(j, r) * Factorial(r));
+  # end;
+
+  num_bijections := function(i, j)
+    local tmp;
+    if j > i then
+      tmp := i;
+      i := j;
+      j := tmp;
+    fi;
+    return cache[i + 1][j + 1];
+  end;
+
+  # Rank so far
+  r := 0;
+  im := [];
+
+  unused_vals := HashSet([1 .. n]);
+
+  for i in [1 .. n] do
+    num_undefined := num_bijections(n - i, n - r);
+    num_each_defined := num_bijections(n - i, n - r - 1);
+    pos := Random(1, num_undefined + (n - r) * num_each_defined);
+    if pos <= num_undefined then
+      im[i] := 0;  # undefined
+    else
+      pos := pos - num_undefined;
+      im[i] := Set(unused_vals)[QuoInt(pos - 1, num_each_defined) + 1];
+      RemoveSet(unused_vals, im[i]);
+      r := r + 1;
+    fi;
+  od;
+  _RANDOM_PPERM_CACHED_VALUES := WeakPointerObj(cache);
+  return PartialPermNC(im);
+end);

gap/semigroups/semipperm.gi

@@ -24,27 +24,33 @@ InstallMethod(SEMIGROUPS_ProcessRandomArgsCons,
 {filt, params} -> SEMIGROUPS_ProcessRandomArgsCons(IsSemigroup, params));
 
 InstallMethod(RandomSemigroupCons, "for IsPartialPermSemigroup and a list",
-[IsPartialPermSemigroup, IsList], {filt, params} ->
-Semigroup(List([1 .. params[1]], i -> RandomPartialPerm(params[2]))));
+[IsPartialPermSemigroup, IsList],
+function(_, params)
+  return Semigroup(List([1 .. params[1]],
+                        i -> UniformRandomPartialPerm(params[2])));
+end);
 
 InstallMethod(RandomMonoidCons, "for IsPartialPermMonoid and a list",
-[IsPartialPermMonoid, IsList], {filt, params} ->
-Monoid(List([1 .. params[1]], i -> RandomPartialPerm(params[2]))));
+[IsPartialPermMonoid, IsList],
+function(_, params)
+  return Monoid(List([1 .. params[1]],
+                     i -> UniformRandomPartialPerm(params[2])));
+end);
 
 InstallMethod(RandomInverseSemigroupCons,
 "for IsPartialPermSemigroup and a list",
 [IsPartialPermSemigroup, IsList],
 function(_, params)
   return InverseSemigroup(List([1 .. params[1]],
-                               i -> RandomPartialPerm(params[2])));
+                               i -> UniformRandomPartialPerm(params[2])));
 end);
 
 InstallMethod(RandomInverseMonoidCons,
 "for IsPartialPermMonoid and a list",
 [IsPartialPermMonoid, IsList],
 function(_, params)
   return InverseMonoid(List([1 .. params[1]],
-                            i -> RandomPartialPerm(params[2])));
+                            i -> UniformRandomPartialPerm(params[2])));
 end);
 
 #############################################################################