Quelle combinat.gi Sprache: unbekannt

#############################################################################
##
## This file is part of GAP, a system for computational discrete algebra.
## This file's authors include Martin Schönert.
##
## Copyright of GAP belongs to its developers, whose names are too numerous
## to list here. Please refer to the COPYRIGHT file for details.
##
## SPDX-License-Identifier: GPL-2.0-or-later
##
## This file contains method for combinatorics.
##

#############################################################################
##
#F Factorial( <n> ) . . . . . . . . . . . . . . . . factorial of an integer
##
InstallGlobalFunction(Factorial, FACTORIAL_INT);

#############################################################################
##
#F GaussianCoefficient( <n>, <k>, <q> ) . . . . . . . . number of subspaces
##
InstallGlobalFunction(GaussianCoefficient,function ( n, k, q )
local gc, i;
 if k < 0 or n<0 or k>n then
 return 0;
 else
 gc:=1;
 for i in [1..k] do
 gc:=gc*(q^(n-i+1)-1)/(q^i-1);
 od;
 return gc;
 fi;
end);

#############################################################################
##
#F Binomial( <n>, <k> ) . . . . . . . . . binomial coefficient of integers
##
InstallGlobalFunction(Binomial, BINOMIAL_INT);

#############################################################################
##
#F Bell( <n> ) . . . . . . . . . . . . . . . . . value of the Bell sequence
##
InstallGlobalFunction(Bell,function ( n )
 local bell, k, i;
 bell := [ 1 ];
 for i in [1..n-1] do
 bell[i+1] := bell[1];
 for k in [0..i-1] do
 bell[i-k] := bell[i-k] + bell[i-k+1];
 od;
 od;
 return bell[1];
end);

#############################################################################
##
#F Stirling1( <n>, <k> ) . . . . . . . . . Stirling number of the first kind
##
InstallGlobalFunction(Stirling1,function ( n, k )
 local sti, i, j;
 if n < k then
 sti := 0;
 elif n = k then
 sti := 1;
 elif n < 0 and k < 0 then
 sti := Stirling2( -k, -n );
 elif k <= 0 then
 sti := 0;
 else
 sti := [ 1 ];
 for j in [2..n-k+1] do
 sti[j] := 0;
 od;
 for i in [1..k] do
 sti[1] := 1;
 for j in [2..n-k+1] do
 sti[j] := (i+j-2) * sti[j-1] + sti[j];
 od;
 od;
 sti := sti[n-k+1];
 fi;
 return sti;
end);

#############################################################################
##
#F Stirling2( <n>, <k> ) . . . . . . . . Stirling number of the second kind
##
## Uses $S_2(n,k) = (-1)^k \sum_{i=1}^{k}{(-1)^i {k \choose i} i^k} / k!$.
##
InstallGlobalFunction(Stirling2,function ( n, k )
 local sti, bin, fib, i;
 if n < k then
 sti := 0;
 elif n = k then
 sti := 1;
 elif n < 0 and k < 0 then
 sti := Stirling1( -k, -n );
 elif k <= 0 then
 sti := 0;
 else
 bin := 1; # (k 0)
 sti := 0; # (-1)^0 (k 0) 0^k
 fib := 1; # 0!
 for i in [1..k] do
 bin := (k-i+1)/i * bin; # (k i) = (k-(i-1))/i (k i-1)
 sti := bin * i^n - sti; # (-1)^i sum (-1)^j (k j) j^k
 fib := fib * i; # i!
 od;
 sti := sti / fib;
 fi;
 return sti;
end);

#############################################################################
##
#F Combinations( <mset> ) . . . . . . set of sorted sublists of a multiset
##
## 'CombinationsA( <mset>, <m>, <n>, <comb>, )' returns the set of all
## combinations of the multiset <mset>, which has size <n>, that begin with
## '<comb>[[1..-1]]'. To do this it finds all elements of <mset> that
## can go at '<comb>[]' and calls itself recursively for each candidate.
## <m>-1 is the position of '<comb>[-1]' in <mset>, so the candidates for
## '<comb>[]' are exactly the elements 'Set( <mset>[[<m>..<n>]] )'.
##
## 'CombinationsK( <mset>, <m>, <n>, <k>, <comb>, )' returns the set of
## all combinations of the multiset <mset>, which has size <n>, that have
## length '+<k>-1', and that begin with '<comb>[[1..-1]]'. To do this
## it finds all elements of <mset> that can go at '<comb>[]' and calls
## itself recursively for each candidate. <m>-1 is the position of
## '<comb>[-1]' in <mset>, so the candidates for '<comb>[]' are
## exactly the elements 'Set( <mset>[<m>..<n>-<k>+1] )'.
##
## 'Combinations' only calls 'CombinationsA' or 'CombinationsK' with initial
## arguments.
##
DeclareGlobalName( "CombinationsA" );
BindGlobal( "CombinationsA", function ( mset, m, n, comb, i )
 local combs, l;
 if m = n+1 then
 comb := ShallowCopy(comb);
 combs := [ comb ];
 else
 comb := ShallowCopy(comb);
 combs := [ ShallowCopy(comb) ];
 for l in [m..n] do
 if l = m or mset[l] <> mset[l-1] then
 comb[i] := mset[l];
 Append( combs, CombinationsA(mset,l+1,n,comb,i+1) );
 fi;
 od;
 fi;
 return combs;
end );

DeclareGlobalName( "CombinationsK" );
BindGlobal( "CombinationsK", function ( mset, m, n, k, comb, i )
 local combs, l;
 if k = 0 then
 comb := ShallowCopy(comb);
 combs := [ comb ];
 else
 combs := [];
 for l in [m..n-k+1] do
 if l = m or mset[l] <> mset[l-1] then
 comb[i] := mset[l];
 Append( combs, CombinationsK(mset,l+1,n,k-1,comb,i+1) );
 fi;
 od;
 fi;
 return combs;
end );

InstallGlobalFunction(Combinations,function ( mset, arg... )
 local combs;
 mset := ShallowCopy(mset); Sort( mset );
 if Length(arg) = 0 then
 combs := CombinationsA( mset, 1, Length(mset), [], 1 );
 elif Length(arg) = 1 then
 combs := CombinationsK( mset, 1, Length(mset), arg[1], [], 1 );
 else
 Error("usage: Combinations( <mset> [, <k>] )");
 fi;
 return combs;
end);

#############################################################################
##
#F IteratorOfCombinations( <mset>[, <k> ] )
#F EnumeratorOfCombinations( <mset> )
##
InstallGlobalFunction(EnumeratorOfCombinations, function(mset)
 local c, max, l, mods, size, els, ElementNumber, NumberElement;
 c := Collected(mset);
 max := List(c, a-> a[2]);
 els := List(c, a-> a[1]);
 l := Length(max);
 mods := max+1;
 size := Product(mods);
 # a combination can contain els[i] from 0 to max[i] times (mods[i]
 # possibilities), we number the combination that contains a[i] times els[i]
 # for all i by n = 1 + sum_i a[i]*m[i] where m[i] = prod_(j<i) mods[i]
 ElementNumber := function(enu, n)
 local comb, res, i, j;
 if n > size then
 Error("Index ", n, " not bound.");
 fi;
 comb := EmptyPlist(l);
 n := n-1;
 for i in [1..l] do
 comb[i] := n mod mods[i];
 n := (n - comb[i])/mods[i];
 od;
 res := [];
 for i in [1..l] do
 for j in [1..comb[i]] do
 Add(res, els[i]);
 od;
 od;
 return res;
 end;
 NumberElement := function(enu, comb)
 local c, d, pos, n, a, i;
 if not IsList(comb) then
 return fail;
 fi;
 c := Collected(comb);
 d := 0*max;
 for a in c do
 pos := PositionSorted(els, a[1]);
 if not IsBound(els[pos]) or els[pos] <> a[1] or a[2] > max[pos] then
 return fail;
 else
 d[pos] := a[2];
 fi;
 od;
 n := 0;
 for i in [l,l-1..1] do
 n := n*mods[i] + d[i];
 od;
 return n+1;
 end;
 return EnumeratorByFunctions(ListsFamily, rec(
 ElementNumber := ElementNumber,
 NumberElement := NumberElement,
 els := els,
 Length := x->size,
 max := max));
end);

BindGlobal("NextIterator_Combinations_set", function(it)
 local res, comb, k, i, len;
 comb := it!.comb;
 if comb = fail then
 Error("No more elements in iterator.");
 fi;
 # first create combination to return
 res := it!.els{comb};
 # now construct indices for next combination
 len := it!.len;
 k := it!.k;
 for i in [1..k] do
 if i = k or comb[i]+1 < comb[i+1] then
 comb[i] := comb[i] + 1;
 comb{[1..i-1]} := [1..i-1];
 break;
 fi;
 od;
 # check if done
 if k = 0 or comb[k] > len then
 it!.comb := fail;
 fi;
 return res;
end);

# helper function to substitute elements described by r!.comb[j],
# j in [1..i] by smallest possible ones
BindGlobal("Distr_Combinations", function(r, i)
 local max, kk, l, comb, j;
 max := r!.max;
 kk := 0;
 l := Length(max);
 comb := r!.comb;
 for j in [1..i] do
 kk := kk + comb[j];
 comb[j] := 0;
 od;
 for i in [1..l] do
 if kk <= max[i] then
 comb[i] := kk;
 break;
 else
 comb[i] := max[i];
 kk := kk - max[i];
 fi;
 od;
end);
BindGlobal("NextIterator_Combinations_mset", function(it)
 local res, comb, l, els, i, j, max;
 if it!.comb = fail then
 Error("No more elements in iterator.");
 fi;
 comb := it!.comb;
 max := it!.max;
 l := Length(comb);
 # first create the combination to return, this is the time critical
 # code which is more efficient in the proper set case above
 res := EmptyPlist(it!.k);
 els := it!.els;
 for i in [1..l] do
 for j in [1..comb[i]] do
 Add(res, els[i]);
 od;
 od;
 # now find next combination if there is one;
 # for this find smallest element which can be substituted by the next
 # larger element and reset the previous ones to the smallest
 # possible ones
 i := 1;
 while i < l and (comb[i] = 0 or comb[i+1] = max[i+1]) do
 i := i+1;
 od;
 if i = l then
 it!.comb := fail;
 else
 comb[i+1] := comb[i+1] + 1;
 comb[i] := comb[i] - 1;
 Distr_Combinations(it, i);
 fi;
 return res;
end);
BindGlobal("IsDoneIterator_Combinations", function(it)
 return it!.comb = fail;
end);
BindGlobal("ShallowCopy_Combinations", function(it)
 return rec(
 NextIterator := it!.NextIterator,
 IsDoneIterator := it!.IsDoneIterator,
 ShallowCopy := it!.ShallowCopy,
 els := it!.els,
 max := it!.max,
 len := it!.len,
 k := it!.k,
 comb := ShallowCopy(it!.comb));
end);
InstallGlobalFunction(IteratorOfCombinations, function(mset, arg...)
 local k, c, max, els, len, comb, NextFunc;
 len := Length(mset);
 if Length(arg) = 0 then
 # case of one argument, call 2-arg version for each k and concatenate
 return ConcatenationIterators(List([0..len], k->
 IteratorOfCombinations(mset, k)));
 fi;
 k := arg[1];
 if k > len then
 return IteratorList([]);
 elif len = 0 then
 return TrivialIterator( [] );
 fi;
 c := Collected(mset);
 max := List(c, a-> a[2]);
 els := List(c, a-> a[1]);
 if Maximum(max) = 1 then
 # in case of a proper set 'mset' we use 'comb' for indices of
 # elements in current combination; this way the generation
 # of the actual combinations is a bit more efficient than below in the
 # general case of a multiset
 comb := [1..k];
 NextFunc := NextIterator_Combinations_set;
 else
 # the general case of a multiset, here 'comb'
 # describes the combination which contains comb[i] times els[i] for all i
 comb := 0*max;
 comb[1] := k;
 # initialize first combination
 Distr_Combinations(rec(comb := comb,max := max),1);
 NextFunc := NextIterator_Combinations_mset;
 fi;
 return IteratorByFunctions(rec(
 NextIterator := NextFunc,
 IsDoneIterator := IsDoneIterator_Combinations,
 ShallowCopy := ShallowCopy_Combinations,
 els := els,
 max := max,
 len := len,
 k := k,
 comb := comb));
end);

#############################################################################
##
#F NrCombinations( <mset> ) . . . . number of sorted sublists of a multiset
##
## 'NrCombinations' just calls 'NrCombinationsSetA', 'NrCombinationsMSetA',
## 'NrCombinationsSetK' or 'NrCombinationsMSetK' depending on the arguments.
##
## 'NrCombinationsSetA' and 'NrCombinationsSetK' use well known identities.
##
## 'NrCombinationsMSetA' and 'NrCombinationsMSetK' call 'NrCombinationsX',
## and return either the sum or the last element of this list.
##
## 'NrCombinationsX' returns the list 'nrs', such that 'nrs[l+1]' is the
## number of combinations of length l. It uses a recursion formula, taking
## more and more of the elements of <mset>.
##
BindGlobal( "NrCombinationsX", function ( mset, k )
 local nrs, nr, cnt, n, l, i;

 # count how often each element appears
 cnt := List( Collected( mset ), pair -> pair[2] );

 # there is one combination of length 0 and no other combination
 # using none of the elements
 nrs := ListWithIdenticalEntries( k+1, 0 );
 nrs[0+1] := 1;

 # take more and more elements
 for n in [1..Length(cnt)] do

 # loop over the possible lengths of combinations
 for l in [k,k-1..0] do

 # compute the number of combinations of length <l>
 # using only the first <n> elements of <mset>
 nr := 0;
 for i in [0..Minimum(cnt[n],l)] do

 # add the number of combinations of length <l>
 # that consist of <l>- of the first <n>-1 elements
 # and copies of the <n>th element
 nr := nr + nrs[l-i+1];

 od;

 nrs[l+1] := nr;

 od;

 od;

 # return the numbers
 return nrs;
end );

BindGlobal( "NrCombinationsSetA", function ( set )
 local nr;
 nr := 2 ^ Size(set);
 return nr;
end );

BindGlobal( "NrCombinationsMSetA", function ( mset )
 local nr;
 nr := Product( Set(mset), i->Number(mset,j->i=j)+1 );
 return nr;
end );

BindGlobal( "NrCombinationsSetK", function ( set, k )
 local nr;
 if k <= Size(set) then
 nr := Binomial( Size(set), k );
 else
 nr := 0;
 fi;
 return nr;
end );

BindGlobal( "NrCombinationsMSetK", function ( mset, k )
 local nr;
 if k <= Length(mset) then
 nr := NrCombinationsX( mset, k )[k+1];
 else
 nr := 0;
 fi;
 return nr;
end );

InstallGlobalFunction(NrCombinations,function ( mset, arg... )
 local nr;
 mset := ShallowCopy(mset); Sort( mset );
 if Length(arg) = 0 then
 if IsSSortedList( mset ) then
 nr := NrCombinationsSetA( mset );
 else
 nr := NrCombinationsMSetA( mset );
 fi;
 elif Length(arg) = 1 then
 if IsSSortedList( mset ) then
 nr := NrCombinationsSetK( mset, arg[1] );
 else
 nr := NrCombinationsMSetK( mset, arg[1] );
 fi;
 else
 Error("usage: NrCombinations( <mset> [, <k>] )");
 fi;
 return nr;
end);

#############################################################################
##
#F Arrangements( <mset> ) . . . . set of ordered combinations of a multiset
##
## 'ArrangementsA( <mset>, <m>, <n>, <comb>, )' returns the set of all
## arrangements of the multiset <mset>, which has size <n>, that begin with
## '<comb>[[1..-1]]'. To do this it finds all elements of <mset> that
## can go at '<comb>[]' and calls itself recursively for each candidate.
## <m> is a boolean list of size <n> that contains 'true' for every element
## of <mset> that we have not yet taken, so the candidates for '<comb>[]'
## are exactly the elements '<mset>[<l>]' such that '<m>[<l>]' is 'true'.
## Some care must be taken to take a candidate only once if it appears more
## than once in <mset>.
##
## 'ArrangementsK( <mset>, <m>, <n>, <k>, <comb>, )' returns the set of
## all arrangements of the multiset <mset>, which has size <n>, that have
## length '+<k>-1', and that begin with '<comb>[[1..-1]]'. To do this
## it finds all elements of <mset> that can go at '<comb>[]' and calls
## itself recursively for each candidate. <m> is a boolean list of size <n>
## that contains 'true' for every element of <mset> that we have not yet
## taken, so the candidates for '<comb>[]' are exactly the elements
## '<mset>[<l>]' such that '<m>[<l>]' is 'true'. Some care must be taken to
## take a candidate only once if it appears more than once in <mset>.
##
## 'Arrangements' only calls 'ArrangementsA' or 'ArrangementsK' with initial
## arguments.
##
DeclareGlobalName( "ArrangementsA" );
BindGlobal( "ArrangementsA", function ( mset, m, n, comb, i )
 local combs, l;
 if i = n+1 then
 comb := ShallowCopy(comb);
 combs := [ comb ];
 else
 comb := ShallowCopy(comb);
 combs := [ ShallowCopy(comb) ];
 for l in [1..n] do
 if m[l] and (l=1 or m[l-1]=false or mset[l]<>mset[l-1]) then
 comb[i] := mset[l];
 m[l] := false;
 Append( combs, ArrangementsA( mset, m, n, comb, i+1 ) );
 m[l] := true;
 fi;
 od;
 fi;
 return combs;
end );

DeclareGlobalName( "ArrangementsK" );
BindGlobal( "ArrangementsK", function ( mset, m, n, k, comb, i )
 local combs, l;
 if k = 0 then
 comb := ShallowCopy(comb);
 combs := [ comb ];
 else
 combs := [];
 for l in [1..n] do
 if m[l] and (l=1 or m[l-1]=false or mset[l]<>mset[l-1]) then
 comb[i] := mset[l];
 m[l] := false;
 Append( combs, ArrangementsK( mset, m, n, k-1, comb, i+1 ) );
 m[l] := true;
 fi;
 od;
 fi;
 return combs;
end );

InstallGlobalFunction(Arrangements,function ( mset, arg... )
 local combs, m;
 mset := ShallowCopy(mset); Sort( mset );
 m := List( mset, ReturnTrue );
 if Length(arg) = 0 then
 combs := ArrangementsA( mset, m, Length(mset), [], 1 );
 elif Length(arg) = 1 then
 combs := ArrangementsK( mset, m, Length(mset), arg[1], [], 1 );
 else
 Error("usage: Arrangements( <mset> [, <k>] )");
 fi;
 return combs;
end);

#############################################################################
##
#F NrArrangements( <mset> ) . number of ordered combinations of a multiset
##
## 'NrArrangements' just calls 'NrArrangementsSetA', 'NrArrangementsMSetA',
## 'NrArrangementsSetK' or 'NrArrangementsMSetK' depending on the arguments.
##
## 'NrArrangementsSetA' and 'NrArrangementsSetK' use well known identities.
##
## 'NrArrangementsMSetA' and 'NrArrangementsMSetK' call 'NrArrangementsX',
## and return either the sum or the last element of this list.
##
## 'NrArrangementsX' returns the list 'nrs', such that 'nrs[l+1]' is the
## number of arrangements of length l. It uses a recursion formula, taking
## more and more of the elements of <mset>.
##
BindGlobal( "NrArrangementsX", function ( mset, k )
 local nrs, nr, cnt, bin, n, l, i;

 # count how often each element appears
 cnt := List( Collected( mset ), pair -> pair[2] );

 # there is one arrangement of length 0 and no other arrangement
 # using none of the elements
 nrs := ListWithIdenticalEntries( k+1, 0 );
 nrs[0+1] := 1;

 # take more and more elements
 for n in [1..Length(cnt)] do

 # loop over the possible lengths of arrangements
 for l in [k,k-1..0] do

 # compute the number of arrangements of length <l>
 # using only the first <n> elements of <mset>
 nr := 0;
 bin := 1;
 for i in [0..Minimum(cnt[n],l)] do

 # add the number of arrangements of length <l>
 # that consist of <l>- of the first <n>-1 elements
 # and copies of the <n>th element
 nr := nr + bin * nrs[l-i+1];
 bin := bin * (l-i) / (i+1);

 od;

 nrs[l+1] := nr;

 od;

 od;

 # return the numbers
 return nrs;
end );

BindGlobal( "NrArrangementsSetA", function ( set )
 local nr, i;
 nr := 0;
 for i in [0..Size(set)] do
 nr := nr + Product([Size(set)-i+1..Size(set)]);
 od;
 return nr;
end );

BindGlobal( "NrArrangementsMSetA", function ( mset, k )
 local nr;
 nr := Sum( NrArrangementsX( mset, k ) );
 return nr;
end );

BindGlobal( "NrArrangementsSetK", function ( set, k )
 local nr;
 if k <= Size(set) then
 nr := Product([Size(set)-k+1..Size(set)]);
 else
 nr := 0;
 fi;
 return nr;
end );

BindGlobal( "NrArrangementsMSetK", function ( mset, k )
 local nr;
 if k <= Length(mset) then
 nr := NrArrangementsX( mset, k )[k+1];
 else
 nr := 0;
 fi;
 return nr;
end );

InstallGlobalFunction(NrArrangements,function ( mset, arg... )
 local nr;
 mset := ShallowCopy(mset); Sort( mset );
 if Length(arg) = 0 then
 if IsSSortedList( mset ) then
 nr := NrArrangementsSetA( mset );
 else
 nr := NrArrangementsMSetA( mset, Length(mset) );
 fi;
 elif Length(arg) = 1 then
 if not (IsInt(arg[1]) and arg[1] >= 0) then
 Error("<k> must be a nonnegative integer");
 fi;
 if IsSSortedList( mset ) then
 nr := NrArrangementsSetK( mset, arg[1] );
 else
 nr := NrArrangementsMSetK( mset, arg[1] );
 fi;
 else
 Error("usage: NrArrangements( <mset> [, <k>] )");
 fi;
 return nr;
end);

#############################################################################
##
#F UnorderedTuples( <set>, <k> ) . . . . set of unordered tuples from a set
##
## 'UnorderedTuplesK( <set>, <n>, <m>, <k>, <tup>, )' returns the set of
## all unordered tuples of the set <set>, which has size <n>, that have
## length '+<k>-1', and that begin with '<tup>[[1..-1]]'. To do this
## it finds all elements of <set> that can go at '<tup>[]' and calls
## itself recursively for each candidate. <m> is the position of
## '<tup>[-1]' in <set>, so the candidates for '<tup>[]' are exactly
## the elements '<set>[[<m>..<n>]]', since we require that unordered tuples
## be sorted.
##
## 'UnorderedTuples' only calls 'UnorderedTuplesK' with initial arguments.
##
DeclareGlobalName( "UnorderedTuplesK" );
BindGlobal( "UnorderedTuplesK", function ( set, n, m, k, tup, i )
 local tups, l;
 if k = 0 then
 tup := ShallowCopy(tup);
 tups := [ tup ];
 else
 tups := [];
 for l in [m..n] do
 tup[i] := set[l];
 Append( tups, UnorderedTuplesK( set, n, l, k-1, tup, i+1 ) );
 od;
 fi;
 return tups;
end );

InstallGlobalFunction(UnorderedTuples,function ( set, k )
 set := Set(set);
 return UnorderedTuplesK( set, Size(set), 1, k, [], 1 );
end);

#############################################################################
##
#F NrUnorderedTuples( <set>, <k> ) . . number unordered of tuples from a set
##
InstallGlobalFunction(NrUnorderedTuples,function ( set, k )
 return Binomial( Size(Set(set))+k-1, k );
end);

#############################################################################
##
#F IteratorOfCartesianProduct( list1, list2, ... )
#F IteratorOfCartesianProduct( list )
##
## All elements of the cartesian product of lists
## <list1>, <list2>, ... are returned in the lexicographic order.
##
BindGlobal( "IsDoneIterator_Cartesian", iter -> ( iter!.next = false ) );

BindGlobal( "NextIterator_Cartesian",
 function( iter )
 local succ, n, sets, res, i, k;
 succ := iter!.next;
 n := iter!.n;
 sets := iter!.sets;
 res := [];
 i := n;
 while i > 0 do
 res[i] := sets[i][succ[i]];
 i := i-1;
 od;

 if succ = iter!.sizes then
 iter!.next := false;
 else
 succ[n] := succ[n] + 1;
 for k in [n,n-1..2] do
 if succ[k] > iter!.sizes[k] then
 succ[k] := 1;
 succ[k-1] := succ[k-1] + 1;
 else
 break;
 fi;
 od;
 fi;

 return res;
 end);

BindGlobal( "ShallowCopy_Cartesian",
 iter -> rec(
 sets := iter!.sets,
 sizes := iter!.sizes,
 n := iter!.n,
 next := ShallowCopy( iter!.next ) ) );

BindGlobal( "IteratorOfCartesianProduct2",
 function( listsets )
 local s, n;
 if not ForAll( listsets, IsListOrCollection ) and ForAll( listsets, IsFinite ) then
 Error( "Each argument must be a finite list or collection" );
 fi;
 if ForAny( listsets, IsEmpty ) then
 return Iterator( [] );
 fi;
 s := List( listsets, Set );
 n := Length( s );
 # from now s is a list of n finite sets
 return IteratorByFunctions(
 rec( IsDoneIterator := IsDoneIterator_Cartesian,
 NextIterator := NextIterator_Cartesian,
 ShallowCopy := ShallowCopy_Cartesian,
 sets := s, # list of sets
 sizes := MakeImmutable( List( s, Size ) ),
 n := n, # number of sets
 next := 0 * [ 1 .. n ] + 1 ) ); # list of 1's
 end);

InstallGlobalFunction( "IteratorOfCartesianProduct",
 function( arg )
 # this mimics usage of functions Cartesian and Cartesian2
 if Length( arg ) = 1 then
 return IteratorOfCartesianProduct2( arg[1] );
 fi;
 return IteratorOfCartesianProduct2( arg );
 end);

BindGlobal( "NumberElement_Cartesian",
function(enum, x)
 local n, mults, colls, sum, pos, i;

 n:=enum!.n;
 mults:=enum!.mults;
 colls:=enum!.colls;

 if Length(x)<>n then
 return fail;
 fi;

 sum:=0;
 for i in [1..n-1] do
 pos:=Position(colls[i], x[i]);
 if pos=fail then
 return fail;
 else
 pos:=pos-1;
 fi;
 sum:=sum+pos*mults[i];
 od;

 pos:=Position(colls[n], x[n]);

 if pos=fail then
 return fail;
 fi;

 return sum+pos;
end);

BindGlobal( "ElementNumber_Cartesian",
function(enum, x)
 local n, mults, out, i, colls;

 if x>Length(enum) then
 return fail;
 fi;

 x:=x-1;

 n:=enum!.n;
 mults:=enum!.mults;
 colls:=enum!.colls;
 out:=EmptyPlist(n);

 for i in [1..n-1] do
 out[i]:=QuoInt(x, mults[i]);
 x:=x-out[i]*mults[i];
 out[i]:=colls[i][out[i]+1];
 od;
 out[n]:=colls[n][x+1];

 return out;
end);

BindGlobal( "EnumeratorOfCartesianProduct2",
 function(colls)
 local new_colls, mults, k, out, i, j;

 if (not ForAll(colls, IsFinite)) or not (ForAll(colls, IsCollection) or
 ForAll(colls, IsEnumeratorByFunctions)) then
 ErrorNoReturn("usage: each argument must be a finite collection or enumerator,");
 fi;

 new_colls:=[];
 for i in [1..Length(colls)] do
 if IsDomain(colls[i]) then
 new_colls[i]:=Enumerator(colls[i]);
 else
 new_colls[i]:=colls[i];
 fi;
 od;

 mults:=List(new_colls, Length);
 for i in [1..Length(new_colls)-1] do
 k:=1;
 for j in [i+1..Length(new_colls)] do
 k:=k*Length(new_colls[j]);
 od;
 mults[i]:=k;
 od;
 mults[Length(new_colls)]:=0;

 out:=EnumeratorByFunctions(ListsFamily,
 rec( NumberElement := NumberElement_Cartesian,
 ElementNumber := ElementNumber_Cartesian,
 mults:=mults,
 n:=Length(colls),
 colls:=new_colls,
 Length:=enum-> Maximum([mults[1],1])*Length(new_colls[1])));
 SetIsFinite(out, true);
 return out;
end);

InstallGlobalFunction( "EnumeratorOfCartesianProduct",
 function( arg )
 # this mimics usage of functions Cartesian and Cartesian2
 if IsEmpty(arg) or ForAny(arg, IsEmpty) then
 return [];
 elif Length( arg ) = 1 then
 return EnumeratorOfCartesianProduct2( arg[1] );
 fi;
 return EnumeratorOfCartesianProduct2( arg );
end);

#############################################################################
##
#F Tuples( <set>, <k> ) . . . . . . . . . set of ordered tuples from a set
##
## 'TuplesK( <set>, <k>, <tup>, )' returns the set of all tuples of the
## set <set> that have length '+<k>-1', and that begin with
## '<tup>[[1..-1]]'. To do this it loops over all elements of <set>,
## puts them at '<tup>[]' and calls itself recursively.
##
## 'Tuples' only calls 'TuplesK' with initial arguments.
##
DeclareGlobalName( "TuplesK" );
BindGlobal( "TuplesK", function ( set, k, tup, i )
 local tups, l;
 if k = 0 then
 tup := ShallowCopy(tup);
 tups := [ tup ];
 else
 tups := [];
 for l in set do
 tup[i] := l;
 Append( tups, TuplesK( set, k-1, tup, i+1 ) );
 od;
 fi;
 return tups;
end );

InstallGlobalFunction(Tuples,function ( set, k )
 set := Set(set);
 return TuplesK( set, k, [], 1 );
end);

#############################################################################
##
#F EnumeratorOfTuples( <set>, <k> )
##
InstallGlobalFunction( EnumeratorOfTuples, function( set, k )
 local enum;

 # Handle some trivial cases first.
 if k = 0 then
 return Immutable( [ [] ] );
 elif IsEmpty( set ) then
 return Immutable( [] );
 fi;

 # Construct the object.
 enum:= EnumeratorByFunctions( CollectionsFamily( FamilyObj( set ) ), rec(
 # Add the functions.
 ElementNumber:= function( enum, n )
 local nn, t, i;
 nn:= n-1;
 t:= [];
 for i in [ 1 .. enum!.k ] do
 t[i]:= RemInt( nn, Length( enum!.set ) ) + 1;
 nn:= QuoInt( nn, Length( enum!.set ) );
 od;
 if nn <> 0 then
 Error( "<enum>[", n, "] must have an assigned value" );
 fi;
 nn:= enum!.set{ Reversed( t ) };
 MakeImmutable( nn );
 return nn;
 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
 return fail;
 fi;
 n:= 0;
 for i in [ 1 .. enum!.k ] do
 n:= Length( enum!.set ) * n + elm[i] - 1;
 od;
 return n+1;
 end,

 Length:= enum -> Length( enum!.set )^enum!.k,

 PrintObj:= function( enum )
 Print( "EnumeratorOfTuples( ", enum!.set, ", ", enum!.k, " )" );
 end,

 # Add the data.
 set:= Set( set ),
 k:= k ) );

 # We know that this enumerator is strictly sorted.
 SetIsSSortedList( enum, true );

 # Return the result.
 return enum;
 end );

#############################################################################
##
#F IteratorOfTuples( <set>, <n> )
##
## All ordered tuples of length <n> of the set <set>
## are returned in lexicographic order.
##
BindGlobal( "IsDoneIterator_Tuples", iter -> ( iter!.next = false ) );

BindGlobal( "NextIterator_Tuples", function( iter )
 local t, m, n, succ, k;

 t := iter!.next;
 m := iter!.m;
 n := iter!.n;

 if t = iter!.last then
 succ := false;
 else
 succ := ShallowCopy( t );
 succ[n] := succ[n] + 1;
 for k in [n,n-1..2] do
 if succ[k] > m then
 succ[k] := succ[k] - m;
 succ[k-1] := succ[k-1] + 1;
 else
 break;
 fi;
 od;
 fi;

 iter!.next:= succ;
 return iter!.set{t};
 end );

BindGlobal( "ShallowCopy_Tuples",
 iter -> rec( m := iter!.m,
 n := iter!.n,
 last := iter!.last,
 set := iter!.set,
 next := ShallowCopy( iter!.next ) ) );

InstallGlobalFunction( "IteratorOfTuples",
 function( s, n )

 if not ( n=0 or IsPosInt( n ) ) then
 Error( "The second argument <n> must be a non-negative integer" );
 fi;

 if not ( IsCollection( s ) and IsFinite( s ) or IsEmpty( s ) and n=0 ) then
 if s = [] then
 return IteratorByFunctions(
 rec( IsDoneIterator := ReturnTrue,
 NextIterator := NextIterator_Tuples,
 ShallowCopy := ShallowCopy_Tuples,
 next := false) );
 else
 Error( "The first argument <s> must be a finite collection or empty" );
 fi;
 fi;
 s := Set(s);
 # from now on s is a finite set and n is its Cartesian power to be enumerated
 return IteratorByFunctions(
 rec( IsDoneIterator := IsDoneIterator_Tuples,
 NextIterator := NextIterator_Tuples,
 ShallowCopy := ShallowCopy_Tuples,
 set := s,
 m := Size(s),
 last := 0 * [1..n] + ~!.m,
 n := n,
 next := 0 * [ 1 .. n ] + 1 ) );
 end );

#############################################################################
##
#F NrTuples( <set>, <k> ) . . . . . . . number of ordered tuples from a set
##
InstallGlobalFunction(NrTuples,function ( set, k )
 return Size(Set(set)) ^ k;
end);

#############################################################################
##
#F PermutationsList( <mset> ) . . . . . . set of permutations of a multiset
##
## 'PermutationsListK( <mset>, <m>, <n>, <k>, <perm>, )' returns the set
## of all permutations of the multiset <mset>, which has size <n>, that
## begin with '<perm>[[1..-1]]'. To do this it finds all elements of
## <mset> that can go at '<perm>[]' and calls itself recursively for each
## candidate. <m> is a boolean list of size <n> that contains 'true' for
## every element of <mset> that we have not yet taken, so the candidates for
## '<perm>[]' are exactly the elements '<mset>[<l>]' such that
## '<m>[<l>]' is 'true'. Some care must be taken to take a candidate only
## once if it appears more than once in <mset>.
##
## 'Permutations' only calls 'PermutationsListK' with initial arguments.
##
DeclareGlobalName( "PermutationsListK" );
BindGlobal( "PermutationsListK", function ( mset, m, n, k, perm, i )
 local perms, l;
 if k = 0 then
 perm := ShallowCopy(perm);
 perms := [ perm ];
 else
 perms := [];
 for l in [1..n] do
 if m[l] and (l=1 or m[l-1]=false or mset[l]<>mset[l-1]) then
 perm[i] := mset[l];
 m[l] := false;
 Append( perms, PermutationsListK(mset,m,n,k-1,perm,i+1) );
 m[l] := true;
 fi;
 od;
 fi;
 return perms;
end );

InstallGlobalFunction(PermutationsList,function ( mset )
 local m;
 mset := ShallowCopy(mset); Sort( mset );
 m := List( mset, ReturnTrue );
 return PermutationsListK(mset,m,Length(mset),Length(mset),[],1);
end);

#############################################################################
##
#F NrPermutationsList( <mset> ) . . . number of permutations of a multiset
##
## 'NrPermutationsList' uses the well known multinomial coefficient formula.
##
InstallGlobalFunction(NrPermutationsList,function ( mset )
 local nr, m;
 nr := Factorial( Length(mset) );
 for m in Set(mset) do
 nr := nr / Factorial( Number( mset, i->i = m ) );
 od;
 return nr;
end);

#############################################################################
##
#F Derangements( <list> ) . . . . set of fixpointfree permutations of a list
##
## 'DerangementsK( <mset>, <m>, <n>, <list>, <k>, <perm>, )' returns the
## set of all permutations of the multiset <mset>, which has size <n>, that
## have no element at the same position as <list>, and that begin with
## '<perm>[[1..-1]]'. To do this it finds all elements of <mset> that
## can go at '<perm>[]' and calls itself recursively for each candidate.
## <m> is a boolean list of size <n> that contains 'true' for every element
## that we have not yet taken, so the candidates for '<perm>[]' are the
## elements '<mset>[<l>]' such that '<m>[<l>]' is 'true'. Some care must be
## taken to take a candidate only once if it append more than once in
## <mset>.
##
DeclareGlobalName( "DerangementsK" );
BindGlobal( "DerangementsK", function ( mset, m, n, list, k, perm, i )
 local perms, l;
 if k = 0 then
 perm := ShallowCopy(perm);
 perms := [ perm ];
 else
 perms := [];
 for l in [1..n] do
 if m[l] and (l=1 or m[l-1]=false or mset[l]<>mset[l-1])
 and mset[l] <> list[i] then
 perm[i] := mset[l];
 m[l] := false;
 Append( perms, DerangementsK(mset,m,n,list,k-1,perm,i+1) );
 m[l] := true;
 fi;
 od;
 fi;
 return perms;
end );

InstallGlobalFunction(Derangements,function ( list )
 local mset, m;
 mset := ShallowCopy(list); Sort( mset );
 m := List( mset, ReturnTrue );
 return DerangementsK(mset,m,Length(mset),list,Length(mset),[],1);
end);

#############################################################################
##
#F NrDerangements( <list> ) . number of fixpointfree permutations of a list
##
## 'NrDerangements' uses well known identities if <mset> is a proper set.
## If <mset> is a multiset it uses 'NrDerangementsK', which works just like
## 'DerangementsK'.
##
DeclareGlobalName( "NrDerangementsK" );
BindGlobal( "NrDerangementsK", function ( mset, m, n, list, k, i )
 local perms, l;
 if k = 0 then
 perms := 1;
 else
 perms := 0;
 for l in [1..n] do
 if m[l] and (l=1 or m[l-1]=false or mset[l]<>mset[l-1])
 and mset[l] <> list[i] then
 m[l] := false;
 perms := perms + NrDerangementsK(mset,m,n,list,k-1,i+1);
 m[l] := true;
 fi;
 od;
 fi;
 return perms;
end );

InstallGlobalFunction(NrDerangements,function ( list )
 local nr, mset, m, i;
 mset := ShallowCopy(list); Sort( mset );
 if IsSSortedList(mset) then
 if Size(mset) = 0 then
 nr := 1;
 elif Size(mset) = 1 then
 nr := 0;
 else
 m := - Factorial(Size(mset));
 nr := 0;
 for i in [2..Size(mset)] do
 m := - m / i;
 nr := nr + m;
 od;
 fi;
 else
 m := List( mset, ReturnTrue );
 nr := NrDerangementsK(mset,m,Length(mset),list,Length(mset),1);
 fi;
 return nr;
end);

#############################################################################
##
#F Permanent( <mat> ) . . . . . . . . . . . . . . . . permanent of a matrix
##
DeclareGlobalName( "Permanent2" );
BindGlobal( "Permanent2", function ( mat, m, n, r, v, i, sum )
 local p, k;
 if i = n+1 then
 p := v;
 for k in sum do p := p * k; od;
 else
 p := Permanent2( mat, m, n, r, v, i+1, sum )
 + Permanent2( mat, m, n, r+1, v*(r-m)/(n-r), i+1, sum+mat[i] );
 fi;
 return p;
end );

InstallMethod(Permanent,
 "for matrices",
 [ IsMatrix ],
function ( mat )
 local m, n;

 m := NrRows(mat);
 n := NrCols(mat);
 while n<m do
 Error("Matrix may not have fewer columns than rows");
 od;
 mat := TransposedMat(mat);
 return Permanent2( mat, m, n, 0, (-1)^m*Binomial(n,m), 1, 0*mat[1] );
end);

#############################################################################
##
#F PartitionsSet( <set> ) . . . . . . . . . . . set of partitions of a set
##
## 'PartitionsSetA( <set>, <n>, <m>, <o>, <part>, , <j> )' returns the
## set of all partitions of the set <set>, which has size <n>, that begin
## with '<part>[[1..-1]]' and where the -th set begins with
## '<part>[][[1..<j>]]'. To do so it does two things. It finds all
## elements of <mset> that can go at '<part>[][<j>+1]' and calls itself
## recursively for each candidate. And it considers the set '<part>[]'
## to be complete and starts a new set '<part>[+1]', which must start
## with the smallest element of <mset> not yet taken, because we require the
## returned partitions to be sorted lexicographically. <mset> is a boolean
## list that contains 'true' for every element of <mset> not yet taken. <o>
## is the position of '<part>[][<j>]' in <mset>, so the candidates for
## '<part>[][<j>+1]' are those elements '<mset>[<l>]' for which '<o> <
## <l>' and '<m>[<l>]' is 'true'.
##
## 'PartitionsSetK( <set>, <n>, <m>, <o>, <k>, <part>, , <j> )' returns
## the set of all partitions of the set <set>, which has size <n>, that have
## '<k>+-1' subsets, and that begin with '<part>[[1..-1]]' and where
## the -th set begins with '<part>[][[1..<j>]]'. To do so it does two
## things. It finds all elements of <mset> that can go at
## '<part>[][<j>+1]' and calls itself recursively for each candidate.
## And, if <k> is larger than 1, it considers the set '<part>[]' to be
## complete and starts a new set '<part>[+1]', which must start with the
## smallest element of <mset> not yet taken, because we require the returned
## partitions to be sorted lexicographically. <mset> is a boolean list that
## contains 'true' for every element of <mset> not yet taken. <o> is the
## position of '<part>[][<j>]' in <mset>, so the candidates for
## '<part>[][<j>+1]' are those elements '<mset>[<l>]' for which '<o> <
## <l>' and '<m>[<l>]' is 'true'.
##
## 'PartitionsSet' only calls 'PartitionsSetA' or 'PartitionsSetK' with
## initial arguments.
##
DeclareGlobalName( "PartitionsSetA" );
BindGlobal( "PartitionsSetA", function ( set, n, m, o, part, i, j )
 local parts, npart, l;
 l := Position(m,true);
 if l = fail then
 part := List(part,ShallowCopy);
 parts := [ part ];
 else
 npart := ShallowCopy(part);
 m[l] := false;
 npart[i+1] := [ set[l] ];
 parts := PartitionsSetA(set,n,m,l+1,npart,i+1,1);
 m[l] := true;
 part := ShallowCopy(part);
 part[i] := ShallowCopy(part[i]);
 for l in [o..n] do
 if m[l] then
 m[l] := false;
 part[i][j+1] := set[l];
 Append( parts, PartitionsSetA(set,n,m,l+1,part,i,j+1));
 m[l] := true;
 fi;
 od;
 fi;
 return parts;
end );

DeclareGlobalName( "PartitionsSetK" );
BindGlobal( "PartitionsSetK", function ( set, n, m, o, k, part, i, j )
 local parts, npart, l;
 l := Position(m,true);
 parts := [];
 if k = 1 then
 part := List(part,ShallowCopy);
 for l in [k..n] do
 if m[l] then
 Add( part[i], set[l] );
 fi;
 od;
 parts := [ part ];
 elif l <> fail then
 npart := ShallowCopy(part);
 m[l] := false;
 npart[i+1] := [ set[l] ];
 parts := PartitionsSetK(set,n,m,l+1,k-1,npart,i+1,1);
 m[l] := true;
 part := ShallowCopy(part);
 part[i] := ShallowCopy(part[i]);
 for l in [o..n] do
 if m[l] then
 m[l] := false;
 part[i][j+1] := set[l];
 Append( parts, PartitionsSetK(set,n,m,l+1,k,part,i,j+1));
 m[l] := true;
 fi;
 od;
 fi;
 return parts;
end );

InstallGlobalFunction(PartitionsSet,function ( set, arg... )
 local parts, m, k;
 if not IsSSortedList(set) then
 Error("PartitionsSet: <set> must be a set");
 fi;
 if Length(arg) = 0 then
 if set = [] then
 parts := [ [ ] ];
 else
 m := List( set, ReturnTrue );
 m[1] := false;
 parts := PartitionsSetA(set,Length(set),m,2,[[set[1]]],1,1);
 fi;
 elif Length(arg) = 1 then
 k := arg[1];
 if set = [] then
 if k = 0 then
 parts := [ [ ] ];
 else
 parts := [ ];
 fi;
 else
 m := List( set, ReturnTrue );
 m[1] := false;
 parts := PartitionsSetK(
 set, Length(set), m, 2, k, [[set[1]]], 1, 1 );
 fi;
 else
 Error("usage: PartitionsSet( <n> [, <k>] )");
 fi;
 return parts;
end);

#############################################################################
##
#F NrPartitionsSet( <set> ) . . . . . . . . . number of partitions of a set
##
InstallGlobalFunction(NrPartitionsSet,function ( set, arg... )
 local nr;
 if not IsSSortedList(set) then
 Error("NrPartitionsSet: <set> must be a set");
 fi;
 if Length(arg) = 0 then
 nr := Bell( Size(set) );
 elif Length(arg) = 1 then
 nr := Stirling2( Size(set), arg[1] );
 else
 Error("usage: NrPartitionsSet( <n> [, <k>] )");
 fi;
 return nr;
end);

#############################################################################
##
#F Partitions( <n> ) . . . . . . . . . . . . set of partitions of an integer
##
## 'PartitionsA( <n>, <m>, <part>, )' returns the set of all partitions
## of '<n> + Sum(<part>[[1..-1]])' that begin with '<part>[[1..-1]]'.
## To do so it finds all values that can go at '<part>[]' and calls
## itself recursively for each candidate. <m> is '<part>[-1]', so the
## candidates for '<part>[]' are '[1..Minimum(<m>,<n>)]', since we
## require that partitions are nonincreasing.
##
## There is one hack that needs some comments. Each call to 'PartitionsA'
## contributes one partition without going into recursion, namely the
## 'Concatenation( <part>[[1..-1]], [1,1,...,1] )'. Of all partitions
## returned by 'PartitionsA' this is the smallest, i.e., it will be the
## first one in the result set. Therefore it is put into the result set
## before anything else is done. However it is not immediately padded with
## 1, this is the last thing 'PartitionsA' does before returning. In the
## meantime the list is used as a temporary that is passed to recursive
## invocations. Note that the fact that each call contributes one partition
## without going into recursion means that the number of recursive calls to
## 'PartitionsA' (and the number of calls to 'ShallowCopy') is equal to
## 'NrPartitions(<n>)'.
##
## 'PartitionsK( <n>, <m>, <k>, <part>, )' returns the set of all
## partitions of '<n> + Sum(<part>[[1..-1]])' that have length
## '<k>+-1' and that begin with '<part>[[1..-1]]'. To do so it finds
## all values that can go at '<part>[]' and calls itself recursively for
## each candidate. <m> is '<part>[-1]', so the candidates for
## '<part>[]' must be less than or equal to <m>, since we require that
## partitions are nonincreasing. Also '<part>[]' must be \<=
## '<n>+1-<k>', since we need at least <k>-1 ones to fill the <k>-1
## positions of <part> remaining after filling '<part>[]'. On the other
## hand '<part>[]' must be >= '<n>/<k>', because otherwise we cannot
## fill the <k>-1 remaining positions nonincreasingly. It is not difficult
## to show that for each candidate satisfying these properties there is
## indeed a partition, i.e., we never run into a dead end.
##
## 'Partitions' only calls 'PartitionsA' or 'PartitionsK' with initial
## arguments.
##
DeclareGlobalName( "PartitionsA" );
BindGlobal( "PartitionsA", function ( n, m, part, i )
 local parts, l;
 if n = 0 then
 part := ShallowCopy(part);
 parts := [ part ];
 elif n <= m then
 part := ShallowCopy(part);
 parts := [ part ];
 for l in [2..n] do
 part[i] := l;
 Append( parts, PartitionsA( n-l, l, part, i+1 ) );
 od;
 for l in [i..i+n-1] do
 part[l] := 1;
 od;
 else
 part := ShallowCopy(part);
 parts := [ part ];
 for l in [2..m] do
 part[i] := l;
 Append( parts, PartitionsA( n-l, l, part, i+1 ) );
 od;
 for l in [i..i+n-1] do
 part[l] := 1;
 od;
 fi;
 return parts;
end );

DeclareGlobalName( "PartitionsK" );
BindGlobal( "PartitionsK", function ( n, m, k, part, i )
 local parts, l;
 if k = 1 then
 part := ShallowCopy(part);
 part[i] := n;
 parts := [ part ];
 elif n+1-k < m then
 parts := [];
 for l in [QuoInt(n+k-1,k)..n+1-k] do
 part[i] := l;
 Append( parts, PartitionsK( n-l, l, k-1, part, i+1 ) );
 od;
 else
 parts := [];
 for l in [QuoInt(n+k-1,k)..m] do
 part[i] := l;
 Append( parts, PartitionsK( n-l, l, k-1, part, i+1 ) );
 od;
 fi;
 return parts;
end );

# The following used to be `Partitions' but was renamed, because
# the new `Partitions' is much faster and produces less garbage, see
# below.
InstallGlobalFunction(PartitionsRecursively,function ( n, arg... )
 local parts, k;
 if Length(arg) = 0 then
 parts := PartitionsA( n, n, [], 1 );
 elif Length(arg) = 1 then
 k := arg[1];
 if n = 0 then
 if k = 0 then
 parts := [ [ ] ];
 else
 parts := [ ];
 fi;
 else
 if k = 0 then
 parts := [ ];
 else
 parts := PartitionsK( n, n, k, [], 1 );
 fi;
 fi;
 else
 Error("usage: Partitions( <n> [, <k>] )");
 fi;
 return parts;
end);

BindGlobal( "GPartitionsEasy", function(n)
 # Returns a list of all Partitions of n, sorted lexicographically.
 # Algorithm/Proof: Let P_n be the set of partitions of n.
 # Let B_n^k be the set of partitions of n with all parts less or equal to k.
 # Then P_n := Union_{k=1}^n [k] + B_{n-k}^k, where "[k]+" means, that
 # a part k is added. Note that the union is a disjoint union.
 # The algorithm first enumerates B_{n-k}^k for k=1,2,...,n-1 and then
 # puts everything together by adding the greatest part.
 # The GAP list B has as its j'th entry B[j] := B_{n-j}^j for j=1,...,n-1.
 # Note the greatest part of all partitions in all of B is less than or
 # equal to QuoInt(n,2).
 # The first stage of the algorithm consists of a loop, where k runs
 # from 1 to QuoInt(n,2) and for each k all partitions are added to all
 # B[j] with greatest part k. Because we run j in descending direction,
 # we already have B[j+k] (partitions of n-j-k) ready up to greatest part k
 # when we handle for B[j] (partitions of n-j) the partitions with greatest
 # part k.
 # In the second stage we only have to add the correct greatest part to get
 # a partition of n.
 # Note that `GPartitions' improves this by including the work for the
 # second step in the first one, such that less garbage objects are generated.
 # n must be a natural number >= 1.
 local B,j,k,l,p,res;
 B := List([1..n-1],x->[]);
 for k in [1..QuoInt(n,2)] do
 # Now we add all partitions for all entries of B with greatest part k.
 Add(B[n-k],[k]); # the trivial partition with greatest part k
 for j in [n-k-1,n-k-2..k] do
 # exactly in those are partitions with greatest part k. Think!
 # we handle B[j] (partitions of n-j) with greatest part k
 for p in B[j+k] do # those are partitions of n-j-k
 l := [k];
 Append(l,p); # This prolonges the bag without creating garbage!
 Add(B[j],l);
 od;
 od;
 od;
 res := []; # here we collect the result
 for k in [1..n-1] do # handle partitions with greatest part k
 for p in B[k] do # use B[k] = B_{n-k}^k
 l := [k]; # add a part k
 Append(l,p);
 Add(res,l); # collect
 od;
 od;
 Add(res,[n]); # one more case
 return res;
end );

BindGlobal( "GPartitions", function(n)
 # Returns a list of all Partitions of n, sorted lexicographically.
 # Algorithm/Proof: See first the comment of `GPartitionsEasy'.
 # This function does exactly the same as `GPartitionsEasy' by the same
 # algorithm, but it produces nearly no garbage, because in contrast
 # to `GPartitionsEasy' the greatest part added in the second stage is
 # already added in the first stage.
 # n must be a natural number >= 1.
 local B,j,k,l,p;
 B := List([1..n],x->[]);
 for k in [1..QuoInt(n,2)] do
 # Now we add all partitions for all entries of B with greatest part k.
 Add(B[n-k],[n-k,k]); # the trivial partition with greatest part k
 for j in [n-k-1,n-k-2..k] do
 # exactly in those are partitions with greatest part k. Think!
 # we handle B[j] (partitions of n-j) with greatest part k
 for p in B[j+k] do # those are partitions of n-j-k
 l := [j]; # This is the greatest part for stage 2
 Append(l,p); # This prolonges the bag without creating garbage!
 l[2] := k; # here used to be the greatest part for stage 2, now k
 Add(B[j],l);
 od;
 od;
 od;
 B[n][1] := [n]; # one more case
 return Concatenation(B);
end );

BindGlobal( "GPartitionsNrPartsHelper", function(n,m,ones)
 # Helper function for GPartitionsNrParts (see below) for the case
 # m > n. This is used only internally if m > QuoInt(n,2), because then
 # the standard routine does not work. Here we just calculate all partitions
 # of n and append a part m to it. We use exactly the algorithm in
 # `GPartitions'.
 local B,j,k,p,res;
 B := List([1..n-1],x->[]);
 for k in [1..QuoInt(n,2)] do
 # Now we add all partitions for all entries of B with greatest part k.
 Add(B[n-k],ones[m]+ones[k]); # the trivial partition with greatest part k
 for j in [n-k-1,n-k-2..k] do
 # exactly in those are partitions with greatest part k. Think!
 # we handle B[j] (partitions of n-j) with greatest part k
 for p in B[j+k] do # those are partitions of n-j-k
 Add(B[j],p + ones[k]);
 od;
 od;
 od;
 res := []; # here we collect the result
 for k in [1..n-1] do # handle partitions with greatest part k
 for p in B[k] do # use B[k] = B_{n-k}^k
 AddRowVector(p,ones[k]);
 Add(res,p); # collect
 od;
 od;
 Add(res,ones[m]+ones[n]); # one more case
 return res;
end );

BindGlobal( "GPartitionsNrParts", function(n,m)
 # This function enumerates the set of all partitions of <n> into exactly
 # <m> parts.
 # We call a partition "admissible", if
 # 0) the sum s of its entries is <= n
 # 1) it has less or equal to m parts
 # 2) let g be its greatest part and k the number of parts,
 # (m-k)*g+s <= n
 # [this means that it may eventually lead to a partition of n with
 # exactly m parts]
 # We proceed in steps. In the first step we write down all admissible
 # partitions with exactly 1 part, sorted by their greatest part.
 # In the t-th step (t from 2 to m-2) we use the partitions from step
 # t-1 to enumerate all admissible partitions with exactly t parts
 # sorted by their greatest part. In step m we add exactly the difference
 # of n and the sum of the entries to get a partition of n.
 #
 # We use the following Lemma: Leaving out the greatest part is a
 # surjective mapping of the set of admissible partitions with k parts
 # to the set of admissible partitions of k-1 parts. Therefore we get
 # every admissible partition with k parts from a partition with k-1
 # parts by adding a part which is greater or equal the greatest part.
 #
 # Note that all our partitions are vectors of length m and until the
 # last step we store n-(the sum) in the first entry.
 #
 local B,BB,i,j,k,p,pos,pp,prototype,t;

 # some special cases:
 if n <= 0 or m < 1 then
 return [];
 elif m = 1 then
 return [[n]];
 fi;
 # from now on we have m >= 2

 prototype := [1..m]*0;

 # Note that there are no admissible partitions of s<n with greatest part
 # greater than QuoInt(n,2) and no one-part-admissible partitions with
 # greatest part greater than QuoInt(n,m):
 # Therefore this is step 1:
 B := [];
 for i in [1..QuoInt(n,m)] do
 B[i] := [ShallowCopy(prototype)];
 B[i][1][1] := n-i; # remember: here is the sum of the parts
 B[i][1][m] := i;
 od;
 for i in [QuoInt(n,m)+1..QuoInt(n,2)] do
 B[i] := [];
 od;

 # Now to steps 2 to m-1:
 for t in [2..m-1] do
 BB := List([1..QuoInt(n,2)],i->[]);
 pos := m+1-t; # here we add a number, this is also number of parts to add
 for j in [1..QuoInt(n,2)] do
 # run through B[j] and add greatest part:
 for p in B[j] do
 # add all possible greatest parts:
 for k in [j+1..QuoInt(p[1],pos)] do
 pp := ShallowCopy(p);
 pp[pos] := k;
 pp[1] := pp[1]-k;
 Add(BB[k],pp);
 od;
 p[pos] := j;
 p[1] := p[1]-j;
 Add(BB[j],p);
 od;
 od;
 B := BB;
 od;

 # In step m we only collect everything (the first entry is already OK!):
 BB := List([1..n-m+1],i->[]);
 for j in [1..Length(B)] do
 for p in B[j] do
 Add(BB[p[1]],p);
 od;
 od;
 return Concatenation(BB);
end );

# The following replaces what is now `PartitionsRecursively':
# It now calls `GPartitions' and friends, which is much faster
# and more environment-friendly because it produces less garbage.
# Thanks to Götz Pfeiffer for the ideas!
InstallGlobalFunction(Partitions,function ( n, arg... )
 local parts, k;
 if not IsInt(n) then
 Error("Partitions: <n> must be an integer");
 fi;
 if Length(arg) = 0 then
 if n <= 0 then
 parts := [[]];
 else
 parts := GPartitions( n );
 fi;
 elif Length(arg) = 1 then
 k := arg[1];
 if not IsInt(k) then
 Error("Partitions: <k> must be an integer");
 fi;
 if n < 0 or k < 0 then
 parts := [];
 else
 if n = 0 then
 if k = 0 then
 parts := [ [ ] ];
 else
 parts := [ ];
 fi;
 else
 if k = 0 then
 parts := [ ];
 else
 parts := GPartitionsNrParts( n, k );
 fi;
 fi;
 fi;
 else
 ErrorNoReturn("usage: Partitions( <n> [, <k>] )");
 fi;
 return parts;
end);

#############################################################################
##
#F NrPartitions( <n> [, <k>] ) . . . . . number of partitions of an integer
##
## To compute $p(n) = NrPartitions(n)$ we use Euler\'s theorem, that asserts
## $p(n) = \sum_{k>0}{ (-1)^{k+1} (p(n-(3m^2-m)/2) + p(n-(3m^2+m)/2)) }$.
##
## To compute $p(n,k)$ we use $p(m,1) = p(m,m) = 1$, $p(m,l) = 0$ if $m\<l$,
## and the recurrence $p(m,l) = p(m-1,l-1) + p(m-l,l)$ if $1 \< l \< m$.
## This recurrence can be proved by splitting the number of ways to write $m$
## as a sum of $l$ summands in two subsets, those sums that have 1 as a
## summand and those that do not. The number of ways to write $m$ as a sum
## of $l$ summands that have 1 as a summand is $p(m-1,l-1)$, because we can
## take away the 1 and obtain a new sums with $l-1$ summands and value
## $m-1$. The number of ways to write $m$ as a sum of $l$ summands such
## that no summand is 1 is $P(m-l,l)$, because we can subtract 1 from each
## summand and obtain new sums that still have $l$ summands but value $m-l$.
##
InstallGlobalFunction(NrPartitions,function ( n, arg... )
 local s, m, p, k, l;

 if Length(arg) = 0 then
 s := 1; # p(0) = 1
 p := [ s ];
 for m in [1..n] do
 s := 0;
 k := 1;
 l := 1; # k*(3*k-1)/2
 while 0 <= m-(l+k) do
 s := s - (-1)^k * (p[m-l+1] + p[m-(l+k)+1]);
 k := k + 1;
 l := l + 3*k - 2;
 od;
 if 0 <= m-l then
 s := s - (-1)^k * p[m-l+1];
 fi;
 p[m+1] := s;
 od;

 elif Length(arg) = 1 then
 k := arg[1];
 if n = k then
 s := 1;
 elif n < k or k = 0 then
 s := 0;
 else
 p := [];
 for m in [1..n] do
 p[m] := 1; # p(m,1) = 1
 od;
 for l in [2..k] do
 for m in [l+1..n-l+1] do
 p[m] := p[m] + p[m-l]; # p(m,l) = p(m,l-1) + p(m-l,l)
 od;
 od;
 s := p[n-k+1];
 fi;

 else
 Error("usage: NrPartitions( <n> [, <k>] )");
 fi;

 return s;
end);

#############################################################################
##
#F PartitionsGreatestLE( <n>, <m> ) . . . set of partitions of n parts <= m
##
## returns the set of all (unordered) partitions of the integer <n> having
## parts less or equal to the integer <m>.
##

BindGlobal( "GPartitionsGreatestLEEasy", function(n,m)
 # Returns a list of all Partitions of n with greatest part less or equal
 # than m, sorted lexicographically.
 # This works essentially as `GPartitions', but the greatest parts are
 # limited.
 # Algorithm/Proof:
 # Let B_n^k be the set of partitions of n with all parts less or equal to k.
 # Then P_n^m := Union_{k=1}^m [k] + B_{n-k}^k}, where "[k]+"
 # means, that a part k is added. Note that the union is a disjoint union.
 # Note that in the end we only need B_{n-k}^k for k<=m but to produce them
 # we need also partial information about B_{n-k}^k for k>m.
 # The algorithm first enumerates B_{n-k}^k for k=1,2,...,m and begins
 # to enumerate B_{n-k}^k for k>m as necessary and then puts everything
 # together by adding the greatest part.
 # The GAP list B has as its j'th entry B[j] := B_{n-j}^j for j=1,...,n-1.
 # Note the greatest part of all partitions in all of B is less than or
 # equal to QuoInt(n,2) and less than or equal to m.
 # The first stage of the algorithm consists of a loop, where k runs
 # from 1 to min(QuoInt(n,2),m) and for each k all partitions are added to all
 # B[j] with greatest part k. Because we run j in descending direction,
 # we already have B[j+k] (partitions of n-j-k) ready up to greatest part k
 # when we handle for B[j] (partitions of n-j) the partitions with greatest
 # part k.
 # In the second stage we only have to add the correct greatest part to get
 # a partition of n.
 # Note that `GPartitionsGreatestLE' improves this by including the
 # work for the second step in the first one, such that less garbage
 # objects are generated.
 # n and m must be a natural numbers >= 1.
 local B,j,k,l,p,res;
 if m >= n then return GPartitions(n); fi; # a special case
 B := List([1..n-1],x->[]);
 for k in [1..Minimum(QuoInt(n,2),m)] do
 # Now we add all partitions for all entries of B with greatest part k.
 Add(B[n-k],[k]); # the trivial partition with greatest part k
 for j in [n-k-1,n-k-2..k] do
 # exactly in those are partitions with greatest part k. Think!
 # we handle B[j] (partitions of n-j) with greatest part k
 for p in B[j+k] do # those are partitions of n-j-k
 l := [k];
 Append(l,p); # This prolonges the bag without creating garbage!
 Add(B[j],l);
 od;
 od;
 od;
 res := []; # here we collect the result
 for k in [1..m] do # handle partitions with greatest part k
 for p in B[k] do # use B[k] = B_{n-k}^k
 l := [k]; # add a part k
 Append(l,p);
 Add(res,l); # collect
 od;
 od;
 return res;
end );

BindGlobal( "GPartitionsGreatestLE", function(n,m)
 # Returns a list of all Partitions of n with greatest part less or equal
 # than m, sorted lexicographically.
 # This works exactly as `GPartitionsGreatestLEEasy', but faster.
 # This is done by doing all the work necessary for step 2 already in step 1.
 # n and m must be a natural numbers >= 1.
 local B,j,k,l,p,res;
 if m >= n then return GPartitions(n); fi; # a special case
 B := List([1..n-1],x->[]);
 for k in [1..Minimum(QuoInt(n,2),m)] do
 # Now we add all partitions for all entries of B with greatest part k.
 Add(B[n-k],[n-k,k]); # the trivial partition with greatest part k
 for j in [n-k-1,n-k-2..k] do
 # exactly in those are partitions with greatest part k. Think!
 # we handle B[j] (partitions of n-j) with greatest part k
 for p in B[j+k] do # those are partitions of n-j-k
 l := [j]; # for step 2
 Append(l,p); # This prolonges the bag without creating garbage!
 l[2] := k; # here we add a new part k
 Add(B[j],l);
 od;
 od;
 od;
 return Concatenation(B{[1..m]});
end );

InstallGlobalFunction( PartitionsGreatestLE,
function(n,m)
 local parts;
 if not(IsInt(n) and IsInt(m)) then
 ErrorNoReturn("usage: PartitionsGreatestLE( <n>, <m> )");
 elif n < 0 or m < 0 then
 parts := [];
 else
 if n = 0 then
 if m = 0 then
 parts := [ [ ] ];
 else
 parts := [ ];
 fi;
 else
 if m = 0 then
 parts := [ ];
 else
 parts := GPartitionsGreatestLE( n, m );
 fi;
 fi;
 fi;
--> --------------------

--> maximum size reached

--> --------------------

Quelle combinat.gi Sprache: unbekannt

[ Dauer der Verarbeitung: 0.28 Sekunden (vorverarbeitet) ]