From 75c6fb2da97340301379e9c4d7a4d818ad5c0dfa Mon Sep 17 00:00:00 2001 From: reiniscirpons Date: Thu, 7 May 2026 07:02:46 +0100 Subject: [PATCH 1/3] First attempt at PartitionsSet enumerstor --- lib/combinat.gd | 32 +++++++-- lib/combinat.gi | 179 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 207 insertions(+), 4 deletions(-) diff --git a/lib/combinat.gd b/lib/combinat.gd index f89bbc550a..cfbde7edbd 100644 --- a/lib/combinat.gd +++ b/lib/combinat.gd @@ -1169,14 +1169,16 @@ DeclareGlobalFunction("NrRestrictedPartitions"); ## DeclareGlobalFunction( "IteratorOfPartitions" ); - ############################################################################# ## -#F IteratorOfPartitionsSet( [, [ ] ] ) +#F IteratorOfPartitionsSet( [, [, ] ] ) +#F EnumeratorOfPartitionsSet( [, [, ] ] ) ## ## <#GAPDoc Label="IteratorOfPartitionsSet"> ## -## +## Iterator and enumerator of unordered set partitions +## +## ## ## ## returns an iterator @@ -1187,12 +1189,34 @@ DeclareGlobalFunction( "IteratorOfPartitions" ); ## then only partitions of size k are computed. ## If k is given and flag is equal to true, ## then only partitions of size at most k are computed. +##

+## returns an enumerator +## (see ) for all unordered partitions of the +## set set. The arguments k and flag function in the +## same manner as for the iterator. +##

+## The ordering of the partitions from these functions can be different and +## also different from the list returned by . +##

+## +## gap> m:=[1..15];; Add(m, 15); +## gap> NrCombinations(m); +## 49152 +## gap> i := 0;; for c in Combinations(m) do i := i+1; od; +## gap> i; +## 49152 +## gap> cm := EnumeratorOfCombinations(m);; +## gap> cm[1000]; +## [ 1, 2, 3, 6, 7, 8, 9, 10 ] +## gap> Position(cm, [1,13,15,15]); +## 36866 +## ## ## ## <#/GAPDoc> ## DeclareGlobalFunction( "IteratorOfPartitionsSet" ); - +DeclareGlobalFunction( "EnumeratorOfPartitionsSet" ); ############################################################################# ## diff --git a/lib/combinat.gi b/lib/combinat.gi index 672992d0d5..f194c72bb8 100644 --- a/lib/combinat.gi +++ b/lib/combinat.gi @@ -2429,6 +2429,185 @@ InstallGlobalFunction( "IteratorOfPartitions", function( n ) next := 0 * [ 1 .. n ] + 1 ) ); end ); +############################################################################# +## +#F EnumeratorOfPartitionsSet( [, [, ] ] ) +## + + +InstallGlobalFunction(EnumeratorOfPartitionsSet , function(s, arg...) + local k, r, l, j, cumulative_stirling, ElementNumber, NumberElement; + + if not IsSet(s) then + Error( "EnumeratorOfPartitionsSet: must be a set" ); + fi; + + l := Length(s); + + # Method for finding n-th k-partition of s, does not + # do argument checking, this is done further down when + # instantiating the enumerators. + ElementNumber := function(enu, n, k) + local res, j, kp, np, t, i; + if l = 0 and n = 0 and k = 0 then + return []; + fi; + + res := List([1 .. k], x -> []); + # kp is the largest index of still active parts. + # Every time kp is decremented, we guarantee that + # res[kp + 1] will not be changed anymore. + kp := k; + # np keeps track of the part of n that we still need + # to process. Use it to avoid recursive calls. + np := n; + for j in [l, l - 1 .. 1] do + # Intuitively, we use identity + # Stirling2(j, k') = Stirling2(j-1, k'-1) + k' * Stirling2(j-1, k') + # to recursively encode partitions in the following manner: + # If n < Stirling2(j-1, k'-1), then we recursively construct + # the (k'-1)-partition corresponding to n and add a new singleton part + # consisting solely of s[j] to it. + # Otherwise, we can write n = Stirling2(j-1, k'-1) + i * Stirling(j-1, k') + n' + # for some numbers i and n' such that + # i < k' and n' < Stirling2(j-1, k'). + # We then recursively construct the k'-partition corresponding to n' and + # add s[j] to the i-th part. + t := Stirling2(j - 1, kp - 1); + if np <= t then + # Add last element to part kp, and + # decrement kp to mark res[kp] as complete. + Add(res[kp], s[j]); + kp := kp - 1; + else + # np and i are the numbers n' and i described above. + np := np - t; + t := Stirling2(j - 1, kp); + i := QuoInt(np - 1, t) + 1; + Add(res[i], s[j]); + np := ((np - 1) mod t) + 1; + fi; + od; + return res; + end; + + NumberElement := function(enu, p, k) + local res, kp, ks, sp, perm, i, j, jp; + + if not IsList(p) or Length(p) <> k then + return fail; + fi; + + if Length(p) = 0 then + if l = 0 and k = 0 then + return 1; + fi; + return fail; + fi; + + # ks[i] is the index corresponding to the + # end of the i-th part of p in the flattened + # version of p (stored as sp) + ks := []; + for kp in [1 .. k] do + if not IsList(p[kp]) then + return fail; + fi; + if kp = 1 then + Add(ks, Length(p[kp])); + else + Add(ks, ks[Length(ks)] + Length(p[kp])); + fi; + od; + sp := Flat(p); + + # Every partition of s must store all of the same elements as s. + perm := PermListList(sp, s); + if perm = fail then + return fail; + fi; + + res := 1; + kp := k; + for j in [l, l - 1 .. 1] do + jp := j^perm; + # i is the part of p that s[j] is in + i := PositionSorted(ks, jp); + + # jp = ks[i] if and only if s[j] is the last element of the i-th + # part of p. According to our encoding convention in ElementNumber, + # this means that the remaining partition is a (kp-1)-partition. + if (jp = ks[i]) then + kp := kp - 1; + else + # Otherwise, we perform the reverse manipulation to the one + # modifying np in ElementNumber + res := res + i * Stirling2(j-1, kp) + Stirling2(j-1, kp-1); + fi; + od; + + return res; + end; + + r := rec(s := s, l := l); + if Length(arg) = 1 or Length(arg) = 2 then + k := arg[1]; + if not IsInt(k) then + Error("EnumeratorOfPartitionsSet: must be an integer"); + fi; + fi; + + if (Length(arg) = 2 and arg[2] = true) or (Length(arg) = 0) then + # k parts or less, also handles all partition case. + if Length(arg) = 0 then + k := l; + else + k := Minimum(k, l); + fi; + + cumulative_stirling := [Stirling2(l, 0)]; + for j in [1 .. k] do + Add(cumulative_stirling, cumulative_stirling[j] + Stirling2(l, j)); + od; + r.cumulative_stirling := cumulative_stirling; + r.ElementNumber := function(enu, n) + local kp; + if n > cumulative_stirling[k + 1] then + Error("Index ", n, " not bound."); + fi; + kp := PositionSorted(cumulative_stirling, n) - 1; + if kp = 0 then + return ElementNumber(enu, n, kp); + fi; + return ElementNumber(enu, n - cumulative_stirling[kp], kp); + end; + r.NumberElement := function(enu, p) + if Length(p) > k then + return fail; + fi; + return NumberElement(enu, p, Length(p)); + end; + r.Length := x -> cumulative_stirling[k + 1]; + elif (Length(arg) = 2 and arg[2] = false) or (Length(arg) = 1) then + r.ElementNumber := function(enu, n) + if n > Stirling2(l, k) then + Error("Index ", n, " not bound."); + fi; + return ElementNumber(enu, n, k); + end; + r.NumberElement := function(enu, p) + return NumberElement(enu, p, k); + end; + r.Length := x -> Stirling2(l, k); + elif Length(arg) = 2 then + Error("EnumeratorOfPartitionsSet: must be true or false"); + else + Error( "usage: EnumeratorOfPartitionsSet( [, [, ] ] )" ); + fi; + + r.k := k; + return EnumeratorByFunctions(ListsFamily, r); +end); ############################################################################# ## From 24897353de7cef00515eee5ea9c443bc3f61a822 Mon Sep 17 00:00:00 2001 From: reiniscirpons Date: Thu, 7 May 2026 15:35:34 +0100 Subject: [PATCH 2/3] Small fix to ordering --- lib/combinat.gi | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/lib/combinat.gi b/lib/combinat.gi index f194c72bb8..0c104b2cdb 100644 --- a/lib/combinat.gi +++ b/lib/combinat.gi @@ -2488,13 +2488,13 @@ InstallGlobalFunction(EnumeratorOfPartitionsSet , function(s, arg...) np := ((np - 1) mod t) + 1; fi; od; - return res; + return List(res, Reversed); end; NumberElement := function(enu, p, k) local res, kp, ks, sp, perm, i, j, jp; - if not IsList(p) or Length(p) <> k then + if not IsSet(p) or Length(p) <> k then return fail; fi; @@ -2509,8 +2509,9 @@ InstallGlobalFunction(EnumeratorOfPartitionsSet , function(s, arg...) # end of the i-th part of p in the flattened # version of p (stored as sp) ks := []; + sp := []; for kp in [1 .. k] do - if not IsList(p[kp]) then + if not IsSet(p[kp]) then return fail; fi; if kp = 1 then @@ -2518,8 +2519,8 @@ InstallGlobalFunction(EnumeratorOfPartitionsSet , function(s, arg...) else Add(ks, ks[Length(ks)] + Length(p[kp])); fi; + Append(sp, p[kp]); od; - sp := Flat(p); # Every partition of s must store all of the same elements as s. perm := PermListList(sp, s); @@ -2534,15 +2535,15 @@ InstallGlobalFunction(EnumeratorOfPartitionsSet , function(s, arg...) # i is the part of p that s[j] is in i := PositionSorted(ks, jp); - # jp = ks[i] if and only if s[j] is the last element of the i-th + # Check passes if and only if s[j] is the first element of the i-th # part of p. According to our encoding convention in ElementNumber, # this means that the remaining partition is a (kp-1)-partition. - if (jp = ks[i]) then + if jp = 1 or (i > 1 and ks[i - 1] + 1 = jp) then kp := kp - 1; else # Otherwise, we perform the reverse manipulation to the one # modifying np in ElementNumber - res := res + i * Stirling2(j-1, kp) + Stirling2(j-1, kp-1); + res := res + (i-1) * Stirling2(j-1, kp) + Stirling2(j-1, kp-1); fi; od; @@ -2585,7 +2586,7 @@ InstallGlobalFunction(EnumeratorOfPartitionsSet , function(s, arg...) if Length(p) > k then return fail; fi; - return NumberElement(enu, p, Length(p)); + return cumulative_stirling[Length(p)] + NumberElement(enu, p, Length(p)); end; r.Length := x -> cumulative_stirling[k + 1]; elif (Length(arg) = 2 and arg[2] = false) or (Length(arg) = 1) then From cc2335427f1e27974c0e526d901cbd06b0217c27 Mon Sep 17 00:00:00 2001 From: reiniscirpons Date: Fri, 8 May 2026 13:20:47 +0100 Subject: [PATCH 3/3] Add tests and iron out small bugs --- lib/combinat.gd | 20 ++--- lib/combinat.gi | 9 ++- tst/testinstall/combinat.tst | 139 ++++++++++++++++++++++++++++++++++- 3 files changed, 156 insertions(+), 12 deletions(-) diff --git a/lib/combinat.gd b/lib/combinat.gd index cfbde7edbd..c583649757 100644 --- a/lib/combinat.gd +++ b/lib/combinat.gd @@ -1199,17 +1199,19 @@ DeclareGlobalFunction( "IteratorOfPartitions" ); ## also different from the list returned by . ##

## -## gap> m:=[1..15];; Add(m, 15); -## gap> NrCombinations(m); -## 49152 -## gap> i := 0;; for c in Combinations(m) do i := i+1; od; +## gap> m:=[1..9];; +## gap> Bell(9); +## 21147 +## gap> i := 0;; for c in PartitionsSet(m) do i := i+1; od; ## gap> i; -## 49152 -## gap> cm := EnumeratorOfCombinations(m);; +## 21147 +## gap> cm := EnumeratorOfPartitionsSet(m);; ## gap> cm[1000]; -## [ 1, 2, 3, 6, 7, 8, 9, 10 ] -## gap> Position(cm, [1,13,15,15]); -## 36866 +## [ [ 1, 2, 4, 9 ], [ 3, 6, 8 ], [ 5, 7 ] ] +## gap> Position(cm, [[1, 3], [2, 5], [4], [6, 7, 8, 9]]); +## 11001 +## gap> cm[11001]; +## [ [ 1, 3 ], [ 2, 5 ], [ 4 ], [ 6, 7, 8, 9 ] ] ## ## ## diff --git a/lib/combinat.gi b/lib/combinat.gi index 0c104b2cdb..0dad5dcfe5 100644 --- a/lib/combinat.gi +++ b/lib/combinat.gi @@ -2583,10 +2583,15 @@ InstallGlobalFunction(EnumeratorOfPartitionsSet , function(s, arg...) return ElementNumber(enu, n - cumulative_stirling[kp], kp); end; r.NumberElement := function(enu, p) - if Length(p) > k then + local t; + if not IsSet(p) or Length(p) > k then return fail; fi; - return cumulative_stirling[Length(p)] + NumberElement(enu, p, Length(p)); + t := NumberElement(enu, p, Length(p)); + if t = fail then + return fail; + fi; + return cumulative_stirling[Length(p)] + t; end; r.Length := x -> cumulative_stirling[k + 1]; elif (Length(arg) = 2 and arg[2] = false) or (Length(arg) = 1) then diff --git a/tst/testinstall/combinat.tst b/tst/testinstall/combinat.tst index f1bf71034f..5b027dac88 100644 --- a/tst/testinstall/combinat.tst +++ b/tst/testinstall/combinat.tst @@ -2,7 +2,7 @@ ## ## This file tests the functions that mainly deal with combinatorics. ## -#@local n,mset,comb1,comb2,comb3,it,pn1,pn2,s,k,x +#@local n,mset,comb1,comb2,comb3,it,pn1,pn2,s,k,x,y gap> START_TEST("combinat.tst"); @@ -549,6 +549,143 @@ gap> for s in [[], [5], [1,2,3,4], [2,5,7], ["a","b","c","d","e"], [3..9]] do > od; > od; +#F EnumeratorOfPartitionsSet( [, [, ] ] ) +gap> EnumeratorOfPartitionsSet(); +Error, Function: number of arguments must be at least 1 (not 0) +gap> EnumeratorOfPartitionsSet(fail); +Error, EnumeratorOfPartitionsSet: must be a set +gap> EnumeratorOfPartitionsSet([],fail); +Error, EnumeratorOfPartitionsSet: must be an integer +gap> EnumeratorOfPartitionsSet([1],1,fail); +Error, EnumeratorOfPartitionsSet: must be true or false +gap> EnumeratorOfPartitionsSet([1],1,true,"too many"); +Error, usage: EnumeratorOfPartitionsSet( [, [, ] ] ) +gap> for s in [ +> [], [5], [1, 2, 3, 4], [2, 5, 7], ["a", "b", "c", "d", "e"], [3 .. 9], +> [[1], [1, 2], [2, 3], [4]], ["a", "ab", "c"], "abcd" +> ] do +> pn1 := PartitionsSet(s); +> pn2 := List(EnumeratorOfPartitionsSet(s)); +> if Length(pn1) <> Length(pn2) then +> Error("wrong number of elements for s = ", s); +> elif Set(pn1) <> Set(pn2) then +> Error("different elements for s = ", s); +> fi; +> for y in pn1 do +> if Position(pn2, y) = fail then +> Error("did not find position of element y = ", y, " when s = ", s); +> fi; +> if pn2[Position(pn2, y)] <> y then +> Error("position does not cancel for y = ", y, " when s = ", s); +> fi; +> od; +> for k in [0 .. Size(s) + 1] do +> pn1 := PartitionsSet(s, k); +> pn2 := List(EnumeratorOfPartitionsSet(s, k)); +> if Length(pn1) <> Length(pn2) then +> Error("wrong number of elements for s = ", s, ", k = ", k); +> elif Set(pn1) <> Set(pn2) then +> Error("different elements for s = ", s, ", k = ", k); +> fi; +> for y in pn1 do +> if Position(pn2, y) = fail then +> Error("did not find position of element y = ", y, " when s = ", s); +> fi; +> if pn2[Position(pn2, y)] <> y then +> Error("position does not cancel for y = ", y, " when s = ", s); +> fi; +> od; +> od; +> for k in [0 .. Size(s) + 1] do +> pn1:= []; +> for x in [0 .. k] do +> Append(pn1, PartitionsSet(s, x)); +> od; +> pn2:= List(IteratorOfPartitionsSet(s, k, true)); +> if Length(pn1) <> Length(pn2) then +> Error("wrong number of elements for s = ", s, ", k <= ", k); +> elif Set(pn1) <> Set(pn2) then +> Error("different elements for s = ", s, ", k <= ", k); +> fi; +> for y in pn1 do +> if Position(pn2, y) = fail then +> Error("did not find position of element y = ", y, " when s = ", s); +> fi; +> if pn2[Position(pn2, y)] <> y then +> Error("position does not cancel for y = ", y, " when s = ", s); +> fi; +> od; +> od; +> od; +gap> pn2 := EnumeratorOfPartitionsSet(["a", "ab", "c", "cd"]);; +gap> Position(pn2, [["a", "ab", "c"], ["cd"]]); +2 +gap> Position(pn2, [["a", "ab"], ["c"], ["cd"]]); +9 +gap> Position(pn2, [["cd"], ["a", "ab", "c"]]); +fail +gap> Position(pn2, [["a", "c", "ab"], ["cd"]]); +fail +gap> Position(pn2, [["a", "ab", "c"], ["d"]]); +fail +gap> Position(pn2, [["a", "ab", "c"], "cd"]); +fail +gap> Position(pn2, [1, 2]); +fail +gap> Position(pn2, 1); +fail +gap> pn2 := EnumeratorOfPartitionsSet(["a", "ab", "c", "cd"], 2);; +gap> Position(pn2, [["a", "ab", "c"], ["cd"]]); +1 +gap> Position(pn2, [["a", "ab"], ["c"], ["cd"]]); +fail +gap> Position(pn2, [["cd"], ["a", "ab", "c"]]); +fail +gap> Position(pn2, [["a", "c", "ab"], ["cd"]]); +fail +gap> Position(pn2, [["a", "ab", "c"], ["d"]]); +fail +gap> Position(pn2, [["a", "ab", "c"], "cd"]); +fail +gap> Position(pn2, [1, 2]); +fail +gap> Position(pn2, 1); +fail +gap> pn2 := EnumeratorOfPartitionsSet(["a", "ab", "c", "cd"], 2, false);; +gap> Position(pn2, [["a", "ab", "c"], ["cd"]]); +1 +gap> Position(pn2, [["a", "ab"], ["c"], ["cd"]]); +fail +gap> Position(pn2, [["cd"], ["a", "ab", "c"]]); +fail +gap> Position(pn2, [["a", "c", "ab"], ["cd"]]); +fail +gap> Position(pn2, [["a", "ab", "c"], ["d"]]); +fail +gap> Position(pn2, [["a", "ab", "c"], "cd"]); +fail +gap> Position(pn2, [1, 2]); +fail +gap> Position(pn2, 1); +fail +gap> pn2 := EnumeratorOfPartitionsSet(["a", "ab", "c", "cd"], 2, true);; +gap> Position(pn2, [["a", "ab", "c"], ["cd"]]); +2 +gap> Position(pn2, [["a", "ab"], ["c"], ["cd"]]); +fail +gap> Position(pn2, [["cd"], ["a", "ab", "c"]]); +fail +gap> Position(pn2, [["a", "c", "ab"], ["cd"]]); +fail +gap> Position(pn2, [["a", "ab", "c"], ["d"]]); +fail +gap> Position(pn2, [["a", "ab", "c"], "cd"]); +fail +gap> Position(pn2, [1, 2]); +fail +gap> Position(pn2, 1); +fail + #F Lucas(

,,) . . . . . . . . . . . . . . value of a lucas sequence gap> Print(List( [0..10], i->Lucas(1,-2,i)[1] ),"\n"); [ 0, 1, 1, 3, 5, 11, 21, 43, 85, 171, 341 ]