Skip to content

Add EnumeratorOfCombinations( mset, k ) and speed up orbits of permutation groups on sets and tuples #6756

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
fingolfin wants to merge 2 commits into
base: master
Choose a base branch
from
+339 −61
Open

Add EnumeratorOfCombinations( mset, k ) and speed up orbits of permutation groups on sets and tuples #6756

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
06b848b
Add k-subset enumerator and use it for orbits on sets
fingolfin Oct 16, 2026
9184750
Check the arguments of EnumeratorOfCombinations
fingolfin Oct 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 11 additions & 9 deletions doc/ref/create.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -492,19 +492,18 @@ function ( set, k )
return nn;
end,
NumberElement := function ( enum, elm )
local n, i;
if not IsList( elm ) then
return fail;
fi;
elm := List( elm, function ( x )
return Position( enum!.set, x );
end );
if fail in elm or Length( elm ) <> enum!.k then
local n, i, p;
if not IsDenseList( elm ) or Length( elm ) <> enum!.k
then
return fail;
fi;
n := 0;
for i in [ 1 .. enum!.k ] do
n := Length( enum!.set ) * n + elm[i] - 1;
p := Position( enum!.set, elm[i] );
if p = fail then
return fail;
fi;
n := Length( enum!.set ) * n + p - 1;
od;
return n + 1;
end,
Expand All @@ -519,6 +518,9 @@ function ( set, k )
set := Set( set ),
k := k ) );
SetIsSSortedList( enum, true );
if IsRange( enum!.set ) or ForAll( enum!.set, IsSmallIntRep ) then
SetFilterObj( enum, IsEnumeratorOfSmallIntLists );
fi;
return enum;
end
]]></Example>
Expand Down
18 changes: 14 additions & 4 deletions lib/combinat.gd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -325,13 +325,13 @@ DeclareGlobalFunction("Combinations");
#############################################################################
##
#F IteratorOfCombinations( mset[, k ] )
#F EnumeratorOfCombinations( mset )
#F EnumeratorOfCombinations( mset[, k ] )
##
## <#GAPDoc Label="IteratorOfCombinations">
## <ManSection>
## <Heading>Iterator and enumerator of combinations</Heading>
## <Func Name="IteratorOfCombinations" Arg='mset[, k]'/>
## <Func Name="EnumeratorOfCombinations" Arg='mset'/>
## <Func Name="EnumeratorOfCombinations" Arg='mset[, k]'/>
##
## <Description>
## <Ref Func="IteratorOfCombinations" /> returns an <Ref Oper="Iterator" />
Expand All @@ -341,8 +341,8 @@ DeclareGlobalFunction("Combinations");
## all combinations.
## <P/>
## <Ref Func="EnumeratorOfCombinations"/> returns an <Ref Attr="Enumerator" />
## of the given multiset <A>mset</A>. Currently only a variant without second
## argument <A>k</A> is implemented.
## for the same combinations.
## With second argument <A>k</A>, it is strictly sorted.
## <P/>
## The ordering of combinations from these functions can be different and also
## different from the list returned by <Ref Func="Combinations"/>.
Expand All @@ -359,6 +359,12 @@ DeclareGlobalFunction("Combinations");
## [ 1, 2, 3, 6, 7, 8, 9, 10 ]
## gap> Position(cm, [1,13,15,15]);
## 36866
## gap> cm := EnumeratorOfCombinations([1..30], 6);
## EnumeratorOfCombinations( [ 1 .. 30 ], 6 )
## gap> Length(cm);
## 593775
## gap> Position(cm, [1,2,3,4,5,7]);
## 2
## </Example>
## </Description>
## </ManSection>
Expand Down Expand Up @@ -639,6 +645,10 @@ DeclareGlobalFunction("Tuples");
##
DeclareGlobalFunction( "EnumeratorOfTuples" );

## Set for the enumerators above if their elements are lists of small
## integers, so that dictionaries can hash these lists by their contents.
DeclareFilter( "IsEnumeratorOfSmallIntLists" );


#############################################################################
##
Expand Down
189 changes: 180 additions & 9 deletions lib/combinat.gi
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -198,10 +198,178 @@ end);
#############################################################################
##
#F IteratorOfCombinations( <mset>[, <k> ] )
#F EnumeratorOfCombinations( <mset> )
#F EnumeratorOfCombinations( <mset>[, <k> ] )
##
InstallGlobalFunction(EnumeratorOfCombinations, function(mset)

# The combinations with <k> entries are numbered in lexicographic order.
# For a set of size n, the combination at positions p_1 < ... < p_k has
# number Binomial(n, k) - Sum_j Binomial(n - p_j, k + 1 - j), using the
# combinatorial number system on the reflected positions n + 1 - p_j.
BindGlobal("ElementNumber_Combinations_set", function(enum, num)
local set, n, r, comb, c, i, lo, hi, mid;
if num > enum!.size then
Error("Index ", num, " not bound, expected at most ", enum!.size);
fi;
set := enum!.set;
n := Length(set);
r := enum!.size - num;
comb := EmptyPlist(enum!.k);
c := n;
for i in [enum!.k, enum!.k-1 .. 1] do
# binary search for the largest lo < c with Binomial(lo, i) <= r
lo := i-1;
hi := c-1;
while lo < hi do
mid := QuoInt(lo+hi+1, 2);
if Binomial(mid, i) <= r then
lo := mid;
else
hi := mid-1;
fi;
od;
Add(comb, set[n-lo]);
r := r - Binomial(lo, i);
c := lo;
od;
return comb;
end);

BindGlobal("NumberElement_Combinations_set", function(enum, comb)
local set, n, k, num, prev, j, p;
k := enum!.k;
if not IsDenseList(comb) or Length(comb) <> k then
return fail;
fi;
set := enum!.set;
n := Length(set);
num := enum!.size;
prev := 0;
for j in [1..k] do
p := Position(set, comb[j]);
if p = fail or p <= prev then
return fail;
fi;
num := num - Binomial(n-p, k+1-j);
prev := p;
od;
return num;
end);

# For a multiset, nr[i][r+1] is the number of combinations with r entries
# from els{[i..Length(els)]}. Of two combinations, the one with more copies
# of the first element in which they differ comes first.
BindGlobal("ElementNumber_Combinations_mset", function(enum, num)
local els, max, nr, r, comb, i, a;
if num > enum!.size then
Error("Index ", num, " not bound, expected at most ", enum!.size);
fi;
els := enum!.els;
max := enum!.max;
nr := enum!.nr;
num := num-1;
r := enum!.k;
comb := [];
for i in [1..Length(els)] do
a := Minimum(max[i], r);
while num >= nr[i+1][r-a+1] do
num := num - nr[i+1][r-a+1];
a := a-1;
od;
Append(comb, ListWithIdenticalEntries(a, els[i]));
r := r-a;
od;
return comb;
end);

BindGlobal("NumberElement_Combinations_mset", function(enum, comb)
local els, max, nr, d, prev, x, p, num, r, i;
if not IsDenseList(comb) or Length(comb) <> enum!.k then
return fail;
fi;
els := enum!.els;
max := enum!.max;
nr := enum!.nr;
d := 0*max;
prev := 1;
for x in comb do
p := Position(els, x);
if p = fail or p < prev or d[p] = max[p] then
return fail;
fi;
d[p] := d[p]+1;
prev := p;
od;
num := 1;
r := enum!.k;
for i in [1..Length(els)] do
num := num + Sum([d[i]+1 .. Minimum(max[i], r)], a -> nr[i+1][r-a+1]);
r := r-d[i];
od;
return num;
end);

BindGlobal("EnumeratorOfCombinations_k", function(mset, k)
local els, c, max, l, nr, i, enum;
if k > Length(mset) then
return Immutable([]);
fi;
els := Set(mset);
if Length(els) = Length(mset) then
enum := EnumeratorByFunctions(CollectionsFamily(FamilyObj(els)), rec(
ElementNumber := ElementNumber_Combinations_set,
NumberElement := NumberElement_Combinations_set,
Length := enum -> enum!.size,
PrintObj := function(enum)
Print("EnumeratorOfCombinations( ");
PrintObj(enum!.set);
Print(", ", enum!.k, " )");
end,
set := els,
k := k,
size := Binomial(Length(els), k)));
else
c := Collected(mset);
max := List(c, a -> a[2]);
l := Length(els);
nr := [];
nr[l+1] := Concatenation([1], ListWithIdenticalEntries(k, 0));
for i in [l, l-1 .. 1] do
nr[i] := List([0..k], r -> Sum([0..Minimum(max[i], r)],
a -> nr[i+1][r-a+1]));
od;
enum := EnumeratorByFunctions(CollectionsFamily(FamilyObj(els)), rec(
ElementNumber := ElementNumber_Combinations_mset,
NumberElement := NumberElement_Combinations_mset,
Length := enum -> enum!.size,
PrintObj := function(enum)
Print("EnumeratorOfCombinations( ");
PrintObj(Concatenation(List([1..Length(enum!.els)],
i -> ListWithIdenticalEntries(enum!.max[i], enum!.els[i]))));
Print(", ", enum!.k, " )");
end,
els := els,
max := max,
nr := nr,
k := k,
size := nr[1][k+1]));
fi;
SetIsSSortedList(enum, true);
if IsRange(els) or ForAll(els, IsSmallIntRep) then
SetFilterObj(enum, IsEnumeratorOfSmallIntLists);
fi;
return enum;
end);

InstallGlobalFunction(EnumeratorOfCombinations, function(mset, arg...)

@james-d-mitchell james-d-mitchell Oct 17, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn't we give an error or something if arg is of length >1?

@james-d-mitchell james-d-mitchell Oct 17, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also shouldn't there be some checks that if arg[1] is present then it is a non-negative integer? At the moment I can do stupid things like:

gap> EnumeratorOfCombinations([1 .. 3], infinity);
[  ]
gap> EnumeratorOfCombinations([1 .. 3], -infinity);
Error, Binomial: <k> must be an integer (not a positional object)
Stack trace:
*[1] Binomial( Length( els ), k )
   @ /Users/jdm3/gap/lib/combinat.gi:329
 [2] EnumeratorOfCombinations_k( mset, arg[1] )
   @ /Users/jdm3/gap/lib/combinat.gi:366
<function "EnumeratorOfCombinations">( <arguments> )
 called from read-eval loop at *stdin*:2
you can enter 'quit;' to quit to outer loop
brk>
gap> EnumeratorOfCombinations([1 .. 3], "a");
[  ]
gap> EnumeratorOfCombinations([1 .. 3], "a", "b");
[  ]

local c, max, l, mods, size, els, ElementNumber, NumberElement;
if Length(arg) > 1 then
Error("usage: EnumeratorOfCombinations( <mset>[, <k>] )");
elif Length(arg) = 1 then
if not (IsInt(arg[1]) and arg[1] >= 0) then
Error("EnumeratorOfCombinations: <k> must be a non-negative integer");
fi;
return EnumeratorOfCombinations_k(mset, arg[1]);
fi;
c := Collected(mset);
max := List(c, a-> a[2]);
els := List(c, a-> a[1]);
Expand Down Expand Up @@ -1006,17 +1174,17 @@ InstallGlobalFunction( EnumeratorOfTuples, function( set, k )
end,

NumberElement:= function( enum, elm )
local n, i;
if not IsList( elm ) then
return fail;
fi;
elm:= List( elm, x -> Position( enum!.set, x ) );
if fail in elm or Length( elm ) <> enum!.k then
local n, i, p;
if not IsDenseList( elm ) or Length( elm ) <> enum!.k then
return fail;
fi;
n:= 0;
for i in [ 1 .. enum!.k ] do
n:= Length( enum!.set ) * n + elm[i] - 1;
p:= Position( enum!.set, elm[i] );
if p = fail then
return fail;
fi;
n:= Length( enum!.set ) * n + p - 1;
od;
return n+1;
end,
Expand All @@ -1033,6 +1201,9 @@ InstallGlobalFunction( EnumeratorOfTuples, function( set, k )

# We know that this enumerator is strictly sorted.
SetIsSSortedList( enum, true );
if IsRange( enum!.set ) or ForAll( enum!.set, IsSmallIntRep ) then
SetFilterObj( enum, IsEnumeratorOfSmallIntLists );
fi;

# Return the result.
return enum;
Expand Down
5 changes: 5 additions & 0 deletions lib/dict.gd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -130,6 +130,11 @@ DeclareCategory( "IsHash", IsLookupDictionary );
##
DeclareGlobalFunction("NewDictionary");

## default for the option blistlimit of NewDictionary: the largest domain
## that gets a position dictionary. Its bit list then takes 4 MB, but the
## value list of a lookup dictionary up to 2^25 words, i.e. 256 MB.
BindGlobal("DICT_BLISTLIMIT", 2^25);

@hulpke hulpke Oct 16, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe indicate how this number arose? Is it 2^17 bytes, i.e. 130kB?


#############################################################################
##
## <#GAPDoc Label="[2]{dict}">
Expand Down
4 changes: 1 addition & 3 deletions lib/dict.gi
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -315,9 +315,7 @@ local hashfun,obj,dom,lookup,maxblist,forcesort;

maxblist:=ValueOption("blistlimit");
if not IsInt(maxblist) then
#2^25 plist (for the lookup list that will hold e.g. positions)
# is 128MB size on a 32-bit system.
maxblist:=2^25;
maxblist:=DICT_BLISTLIMIT;
fi;

# are we given a domain, which can index very quickly?
Expand Down
34 changes: 6 additions & 28 deletions lib/dicthf.gi
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,41 +148,19 @@ end);

#############################################################################
##
#M SparseIntKey(<objcol>)
#M SparseIntKey( <enum>, <list> )
##
InstallMethod(SparseIntKey,"for bounded tuples",true,
[ IsList,IsList and IsCyclotomicCollection ], 0,
InstallMethod(SparseIntKey,"for enumerators of lists of small integers",true,
[ IsEnumeratorOfSmallIntLists, IsObject ], 0,
function(m, v)
if Length(m)<> 3 or m[1]<>"BoundedTuples" then
TryNextMethod();
fi;
# Due to the way BoundedTuples are presently implemented we expect the input
# to the hash function to always be a list of positive immediate integers. This means
# that using HashKeyWholeBag should be safe.
# the entries are immediate, so equal lists have equal entry slots
return function(x)
Assert(1, IsPositionsList(x));
if not IsPlistRep(x) then
x := PlainListCopy(x);
fi;
return HashKeyBag(x, 1, 0, Length(x) * GAPInfo.BytesPerVariable);
return HashKeyBag(x, 1, GAPInfo.BytesPerVariable,
Length(x) * GAPInfo.BytesPerVariable);
end;

# alternative code w/o HashKeyBag
## build a weight vector to distinguish lists. Make entries large while staying clearly within
## immediate int (2^55 replacing 2^60, since we take subsequent primes).
#step:=NextPrimeInt(QuoInt(2^55,Maximum(m[2])*m[3]));
#weights:=[1];
#len:=Length(v);
## up to 56 full, then increasingly reduce
#len:=Minimum(len,8*RootInt(len));
#while Length(weights)<len do
# Add(weights,Last(weights)+step);
# step:=NextPrimeInt(step);
#od;
#return function(a)
# return a*weights;
#end;

end);

BindGlobal( "SparseIntKeyVecListAndMatrix", function(d,m)
Expand Down
Loading
Loading