Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
| Download
GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
Project: cocalc-sagemath-dev-slelievre
Views: 418346#############################################################################
##
#W combinat.gi GAP library Martin Schönert
##
##
#Y Copyright (C) 1996, Lehrstuhl D für Mathematik, RWTH Aachen, Germany
#Y (C) 1998 School Math and Comp. Sci., University of St Andrews, Scotland
#Y Copyright (C) 2002 The GAP Group
##
## This file contains method for combinatorics.
##
#############################################################################
##
#F Factorial( <n> ) . . . . . . . . . . . . . . . . factorial of an integer
##
# can be much further improved, together with Binomial ... (FL)
# but for the moment this is huge improvement over Product([1..n]) for large n
# Factorial(1000000) is no problem now
InstallGlobalFunction(Factorial,function ( n )
local pr;
if n < 0 then Error("<n> must be nonnegative"); fi;
pr := function(l, i, j)
local bound, len, res, l2, k;
bound := 30;
len := j+1-i;
if len < bound then
res := 1;
for k in [i..j] do
res := res*l[k];
od;
return res;
fi;
l2 := QuoInt(len,2);
return pr(l,i,i+l2)*pr(l,i+l2+1,j);
end;
return pr( [1..n], 1, n );
end);
#############################################################################
##
#F Binomial( <n>, <k> ) . . . . . . . . . binomial coefficient of integers
##
InstallGlobalFunction(GaussianCoefficient,function ( n, k, q )
local gc, i, j;
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,function ( n, k )
local bin, i, j;
if k < 0 then
bin := 0;
elif k = 0 then
bin := 1;
elif n < 0 then
bin := (-1)^k * Binomial( -n+k-1, k );
elif n < k then
bin := 0;
elif n = k then
bin := 1;
elif n-k < k then
bin := Binomial( n, n-k );
else
bin := 1; j := 1;
# note that all intermediate results are binomial coefficients itself
# hence integers!
# slight improvement by Frank and Max.
for i in [0..k-1] do
bin := bin * (n-i) / j;
j := j + 1;
od;
fi;
return bin;
end);
#############################################################################
##
#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>, <i> )' returns the set of all
## combinations of the multiset <mset>, which has size <n>, that begin with
## '<comb>[[1..<i>-1]]'. To do this it finds all elements of <mset> that
## can go at '<comb>[<i>]' and calls itself recursively for each candidate.
## <m>-1 is the position of '<comb>[<i>-1]' in <mset>, so the candidates for
## '<comb>[<i>]' are exactly the elements 'Set( <mset>[[<m>..<n>]] )'.
##
## 'CombinationsK( <mset>, <m>, <n>, <k>, <comb>, <i> )' returns the set of
## all combinations of the multiset <mset>, which has size <n>, that have
## length '<i>+<k>-1', and that begin with '<comb>[[1..<i>-1]]'. To do this
## it finds all elements of <mset> that can go at '<comb>[<i>]' and calls
## itself recursively for each candidate. <m>-1 is the position of
## '<comb>[<i>-1]' in <mset>, so the candidates for '<comb>[<i>]' are
## exactly the elements 'Set( <mset>[<m>..<n>-<k>+1] )'.
##
## 'Combinations' only calls 'CombinationsA' or 'CombinationsK' with initial
## arguments.
##
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;
MakeReadOnlyGlobal( "CombinationsA" );
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;
MakeReadOnlyGlobal( "CombinationsK" );
InstallGlobalFunction(Combinations,function ( arg )
local combs, mset;
if Length(arg) = 1 then
mset := ShallowCopy(arg[1]); Sort( mset );
combs := CombinationsA( mset, 1, Length(mset), [], 1 );
elif Length(arg) = 2 then
mset := ShallowCopy(arg[1]); Sort( mset );
combs := CombinationsK( mset, 1, Length(mset), arg[2], [], 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(arg)
local mset, k, c, max, els, len, comb, NextFunc;
mset := arg[1];
len := Length(mset);
if Length(arg) = 1 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[2];
if k > Length(mset) then
return IteratorList([]);
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>-<i> of the first <n>-1 elements
# and <i> 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, k )
local nr;
nr := 2 ^ Size(set);
return nr;
end );
BindGlobal( "NrCombinationsMSetA", function ( mset, k )
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 ( arg )
local nr, mset;
if Length(arg) = 1 then
mset := ShallowCopy(arg[1]); Sort( mset );
if IsSSortedList( mset ) then
nr := NrCombinationsSetA( mset, Length(mset) );
else
nr := NrCombinationsMSetA( mset, Length(mset) );
fi;
elif Length(arg) = 2 then
mset := ShallowCopy(arg[1]); Sort( mset );
if IsSSortedList( mset ) then
nr := NrCombinationsSetK( mset, arg[2] );
else
nr := NrCombinationsMSetK( mset, arg[2] );
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>, <i> )' returns the set of all
## arrangements of the multiset <mset>, which has size <n>, that begin with
## '<comb>[[1..<i>-1]]'. To do this it finds all elements of <mset> that
## can go at '<comb>[<i>]' 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>[<i>]'
## 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>, <i> )' returns the set of
## all arrangements of the multiset <mset>, which has size <n>, that have
## length '<i>+<k>-1', and that begin with '<comb>[[1..<i>-1]]'. To do this
## it finds all elements of <mset> that can go at '<comb>[<i>]' 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>[<i>]' 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.
##
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;
MakeReadOnlyGlobal( "ArrangementsA" );
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;
MakeReadOnlyGlobal( "ArrangementsK" );
InstallGlobalFunction(Arrangements,function ( arg )
local combs, mset, m;
if Length(arg) = 1 then
mset := ShallowCopy(arg[1]); Sort( mset );
m := List( mset, i->true );
combs := ArrangementsA( mset, m, Length(mset), [], 1 );
elif Length(arg) = 2 then
mset := ShallowCopy(arg[1]); Sort( mset );
m := List( mset, i->true );
combs := ArrangementsK( mset, m, Length(mset), arg[2], [], 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>-<i> of the first <n>-1 elements
# and <i> 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, k )
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 ( arg )
local nr, mset;
if Length(arg) = 1 then
mset := ShallowCopy(arg[1]); Sort( mset );
if IsSSortedList( mset ) then
nr := NrArrangementsSetA( mset, Length(mset) );
else
nr := NrArrangementsMSetA( mset, Length(mset) );
fi;
elif Length(arg) = 2 then
if not (IsInt(arg[2]) and arg[2] >= 0) then
Error("<k> must be a nonnegative integer");
fi;
mset := ShallowCopy(arg[1]); Sort( mset );
if IsSSortedList( mset ) then
nr := NrArrangementsSetK( mset, arg[2] );
else
nr := NrArrangementsMSetK( mset, arg[2] );
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>, <i> )' returns the set of
## all unordered tuples of the set <set>, which has size <n>, that have
## length '<i>+<k>-1', and that begin with '<tup>[[1..<i>-1]]'. To do this
## it finds all elements of <set> that can go at '<tup>[<i>]' and calls
## itself recursively for each candidate. <m> is the position of
## '<tup>[<i>-1]' in <set>, so the candidates for '<tup>[<i>]' are exactly
## the elements '<set>[[<m>..<n>]]', since we require that unordered tuples
## be sorted.
##
## 'UnorderedTuples' only calls 'UnorderedTuplesK' with initial arguments.
##
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;
MakeReadOnlyGlobal( "UnorderedTuplesK" );
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(
sizes := iter!.sizes,
n := iter!.n,
next := ShallowCopy( iter!.next ) ) );
BindGlobal( "IteratorOfCartesianProduct2",
function( listsets )
local s, n, x;
if not ForAll( listsets, IsCollection ) and ForAll( listsets, IsFinite ) then
Error( "Each arguments must be a finite collection" );
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 := List( s, Size ), # sizes of sets
n := n, # number of sets
nextelts := List( s, x -> x[1] ), # list of 1st elements
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] );
else
return IteratorOfCartesianProduct2( arg );
fi;
return;
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
Error("usage: each argument must be a finite collection or enumerator,");
return;
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 EmptyPlist(0);
elif Length( arg ) = 1 then
return EnumeratorOfCartesianProduct2( arg[1] );
else
return EnumeratorOfCartesianProduct2( arg );
fi;
return;
end);
#############################################################################
##
#F Tuples( <set>, <k> ) . . . . . . . . . set of ordered tuples from a set
##
## 'TuplesK( <set>, <k>, <tup>, <i> )' returns the set of all tuples of the
## set <set> that have length '<i>+<k>-1', and that begin with
## '<tup>[[1..<i>-1]]'. To do this it loops over all elements of <set>,
## puts them at '<tup>[<i>]' and calls itself recursively.
##
## 'Tuples' only calls 'TuplesK' with initial arguments.
##
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;
MakeReadOnlyGlobal( "TuplesK" );
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>, <i> )' returns the set
## of all permutations of the multiset <mset>, which has size <n>, that
## begin with '<perm>[[1..<i>-1]]'. To do this it finds all elements of
## <mset> that can go at '<perm>[<i>]' 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>[<i>]' 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 apears more than once in <mset>.
##
## 'Permutations' only calls 'PermutationsListK' with initial arguments.
##
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;
MakeReadOnlyGlobal( "PermutationsListK" );
InstallGlobalFunction(PermutationsList,function ( mset )
local m;
mset := ShallowCopy(mset); Sort( mset );
m := List( mset, i->true );
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>, <i> )' 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..<i>-1]]'. To do this it finds all elements of <mset> that
## can go at '<perm>[<i>]' 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>[<i>]' 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>.
##
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;
MakeReadOnlyGlobal( "DerangementsK" );
InstallGlobalFunction(Derangements,function ( list )
local mset, m;
mset := ShallowCopy(list); Sort( mset );
m := List( mset, i->true );
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'.
##
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;
MakeReadOnlyGlobal( "NrDerangementsK" );
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, i->true );
nr := NrDerangementsK(mset,m,Length(mset),list,Length(mset),1);
fi;
return nr;
end);
#############################################################################
##
#F Permanent( <mat> ) . . . . . . . . . . . . . . . . permanent of a matrix
##
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;
MakeReadOnlyGlobal( "Permanent2" );
InstallGlobalFunction(Permanent,function ( mat )
local m, n;
m := Length(mat);
n := Length(mat[1]);
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>, <i>, <j> )' returns the
## set of all partitions of the set <set>, which has size <n>, that begin
## with '<part>[[1..<i>-1]]' and where the <i>-th set begins with
## '<part>[<i>][[1..<j>]]'. To do so it does two things. It finds all
## elements of <mset> that can go at '<part>[<i>][<j>+1]' and calls itself
## recursively for each candidate. And it considers the set '<part>[<i>]'
## to be complete and starts a new set '<part>[<i>+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>[<i>][<j>]' in <mset>, so the candidates for
## '<part>[<i>][<j>+1]' are those elements '<mset>[<l>]' for which '<o> <
## <l>' and '<m>[<l>]' is 'true'.
##
## 'PartitionsSetK( <set>, <n>, <m>, <o>, <k>, <part>, <i>, <j> )' returns
## the set of all partitions of the set <set>, which has size <n>, that have
## '<k>+<i>-1' subsets, and that begin with '<part>[[1..<i>-1]]' and where
## the <i>-th set begins with '<part>[<i>][[1..<j>]]'. To do so it does two
## things. It finds all elements of <mset> that can go at
## '<part>[<i>][<j>+1]' and calls itself recursively for each candidate.
## And, if <k> is larger than 1, it considers the set '<part>[<i>]' to be
## complete and starts a new set '<part>[<i>+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>[<i>][<j>]' in <mset>, so the candidates for
## '<part>[<i>][<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.
##
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;
MakeReadOnlyGlobal( "PartitionsSetA" );
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;
MakeReadOnlyGlobal( "PartitionsSetK" );
InstallGlobalFunction(PartitionsSet,function ( arg )
local parts, set, m;
if Length(arg) = 1 then
set := arg[1];
if not IsSSortedList(arg[1]) then
Error("PartitionsSet: <set> must be a set");
fi;
if set = [] then
parts := [ [ ] ];
else
m := List( set, i->true );
m[1] := false;
parts := PartitionsSetA(set,Length(set),m,2,[[set[1]]],1,1);
fi;
elif Length(arg) = 2 then
set := arg[1];
if not IsSSortedList(set) then
Error("PartitionsSet: <set> must be a set");
fi;
if set = [] then
if arg[2] = 0 then
parts := [ [ ] ];
else
parts := [ ];
fi;
else
m := List( set, i->true );
m[1] := false;
parts := PartitionsSetK(
set, Length(set), m, 2, arg[2], [[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 ( arg )
local nr, set;
if Length(arg) = 1 then
set := arg[1];
if not IsSSortedList(arg[1]) then
Error("NrPartitionsSet: <set> must be a set");
fi;
nr := Bell( Size(set) );
elif Length(arg) = 2 then
set := arg[1];
if not IsSSortedList(set) then
Error("NrPartitionsSet: <set> must be a set");
fi;
nr := Stirling2( Size(set), arg[2] );
else
Error("usage: NrPartitionsSet( <n> [, <k>] )");
fi;
return nr;
end);
#############################################################################
##
#F Partitions( <n> ) . . . . . . . . . . . . set of partitions of an integer
##
## 'PartitionsA( <n>, <m>, <part>, <i> )' returns the set of all partitions
## of '<n> + Sum(<part>[[1..<i>-1]])' that begin with '<part>[[1..<i>-1]]'.
## To do so it finds all values that can go at '<part>[<i>]' and calls
## itself recursively for each candidate. <m> is '<part>[<i>-1]', so the
## candidates for '<part>[<i>]' 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..<i>-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. Therefor 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 befor 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>, <i> )' returns the set of all
## partitions of '<n> + Sum(<part>[[1..<i>-1]])' that have length
## '<k>+<i>-1' and that begin with '<part>[[1..<i>-1]]'. To do so it finds
## all values that can go at '<part>[<i>]' and calls itself recursively for
## each candidate. <m> is '<part>[<i>-1]', so the candidates for
## '<part>[<i>]' must be less than or equal to <m>, since we require that
## partitions are nonincreasing. Also '<part>[<i>]' 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>[<i>]'. On the other
## hand '<part>[<i>]' must be >= '<n>/<k>', because otherwise we can not
## 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.
##
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;
MakeReadOnlyGlobal( "PartitionsA" );
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;
MakeReadOnlyGlobal( "PartitionsK" );
# 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 ( arg )
local parts;
if Length(arg) = 1 then
parts := PartitionsA( arg[1], arg[1], [], 1 );
elif Length(arg) = 2 then
if arg[1] = 0 then
if arg[2] = 0 then
parts := [ [ ] ];
else
parts := [ ];
fi;
else
if arg[2] = 0 then
parts := [ ];
else
parts := PartitionsK( arg[1], arg[1], arg[2], [], 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 ( arg )
local parts;
if Length(arg) = 1 then
if not(IsInt(arg[1])) then
Error("usage: Partitions( <n> [, <k>] )");
else
if arg[1] <= 0 then
parts := [[]];
else
parts := GPartitions( arg[1] );
fi;
fi;
elif Length(arg) = 2 then
if not(IsInt(arg[1]) and IsInt(arg[2])) then
Error("usage: Partitions( <n> [, <k>] )");
return;
elif arg[1] < 0 or arg[2] < 0 then
parts := [];
else
if arg[1] = 0 then
if arg[2] = 0 then
parts := [ [ ] ];
else
parts := [ ];
fi;
else
if arg[2] = 0 then
parts := [ ];
else
parts := GPartitionsNrParts( arg[1], arg[2] );
fi;
fi;
fi;
else
Error("usage: Partitions( <n> [, <k>] )");
return;
fi;
return parts;
end);
#############################################################################
##
#F NrPartitions( <n> ) . . . . . . . . . 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 spliting 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 ( arg )
local s, n, m, p, k, l;
if Length(arg) = 1 then
n := arg[1];
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) = 2 then
if arg[1] = arg[2] then
s := 1;
elif arg[1] < arg[2] or arg[2] = 0 then
s := 0;
else
n := arg[1]; k := arg[2];
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
Error("usage: PartitionsGreatestLE( <n>, <m> )");
return;
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;
return parts;
end);
#############################################################################
##
#F PartitionsGreatestEQ( <n>, <m> ) . . . . set of partitions of n parts = n
##
## returns the set of all (unordered) partitions of the integer <n> having
## greatest part equal to the integer <m>.
##
BindGlobal( "GPartitionsGreatestEQHelper", function(n,m)
# Helper function for GPartitionsGreatestEQ (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,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],[m,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 := [m]; # the greatest part
Append(l,p); # This prolonges the bag without creating garbage!
l[2] := j; # This is the greatest part for stage 2
l[3] := k; # here used to be the greatest part for stage 2, now k
Add(B[j],l);
od;
od;
od;
B[n][1] := [m,n]; # one more case
return Concatenation(B);
end );
BindGlobal( "GPartitionsGreatestEQ", function(n,m)
# Returns a list of all Partitions of n with greatest part equal to
# m, sorted lexicographically.
# This works exactly as `GPartitionsGreatestLE' for n-m and m and
# adds a part m to all partitions. This is however done effectively
# during the work.
# This is the same as `Partitions(n,m)' in the GAP library.
# n and m must be a natural numbers >= 1.
local B,j,k,l,p,res;
if m > n then return []; fi; # a special case
if m = n then return [[m]]; fi; # another special case
n := n - m; # this is >= 1
if m >= n then return GPartitionsGreatestEQHelper(n,m); fi;
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],[m,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 := [m]; # the greatest part m
Append(l,p); # This prolonges the bag without creating garbage!
l[2] := j; # for step 2
l[3] := k; # here we add a new part k
Add(B[j],l);
od;
od;
od;
return Concatenation(B{[1..m]});
end );
InstallGlobalFunction( PartitionsGreatestEQ,
function(n,m)
local parts;
if not(IsInt(n) and IsInt(m)) then
Error("usage: PartitionsGreatestEQ( <n>, <m> )");
return;
elif n < 0 or m < 0 then
parts := [];
else
if m = 0 or n = 0 then
parts := [];
else
parts := GPartitionsGreatestEQ( n, m );
fi;
fi;
return parts;
end);
#############################################################################
##
#F OrderedPartitions( <n> ) . . . . set of ordered partitions of an integer
##
## 'OrderedPartitionsA( <n>, <part>, <i> )' returns the set of all ordered
## partitions of '<n> + Sum(<part>[[1..<i>-1]])' that begin with
## '<part>[[1..<i>-1]]'. To do so it puts all possible values at
## '<part>[<i>]', which are of course exactly the elements of '[1..n]', and
## calls itself recursively.
##
## 'OrderedPartitionsK( <n>, <k>, <part>, <i> )' returns the set of all
## ordered partitions of '<n> + Sum(<part>[[1..<i>-1]])' that have length
## '<k>+<i>-1', and that begin with '<part>[[1..<i>-1]]'. To do so it puts
## all possible values at '<part>[<i>]', which are of course exactly the
## elements of '[1..<n>-<k>+1]', and calls itself recursively.
##
## 'OrderedPartitions' only calls 'OrderedPartitionsA' or
## 'OrderedPartitionsK' with initial arguments.
##
OrderedPartitionsA := function ( n, part, i )
local parts, l;
if n = 0 then
part := ShallowCopy(part);
parts := [ part ];
else
part := ShallowCopy(part);
parts := [];
for l in [1..n-1] do
part[i] := l;
Append( parts, OrderedPartitionsA( n-l, part, i+1 ) );
od;
part[i] := n;
Add( parts, part );
fi;
return parts;
end;
MakeReadOnlyGlobal( "OrderedPartitionsA" );
OrderedPartitionsK := function ( n, k, part, i )
local parts, l;
if k = 1 then
part := ShallowCopy(part);
part[i] := n;
parts := [ part ];
else
parts := [];
for l in [1..n-k+1] do
part[i] := l;
Append( parts, OrderedPartitionsK( n-l, k-1, part, i+1 ) );
od;
fi;
return parts;
end;
MakeReadOnlyGlobal( "OrderedPartitionsK" );
InstallGlobalFunction(OrderedPartitions,function ( arg )
local parts;
if Length(arg) = 1 then
parts := OrderedPartitionsA( arg[1], [], 1 );
elif Length(arg) = 2 then
if arg[1] = 0 then
if arg[2] = 0 then
parts := [ [ ] ];
else
parts := [ ];
fi;
else
if arg[2] = 0 then
parts := [ ];
else
parts := OrderedPartitionsK( arg[1], arg[2], [], 1 );
fi;
fi;
else
Error("usage: OrderedPartitions( <n> [, <k>] )");
fi;
return parts;
end);
#############################################################################
##
#F NrOrderedPartitions( <n> ) . . number of ordered partitions of an integer
##
## 'NrOrderedPartitions' uses well known identities to compute the number of
## ordered partitions of <n>.
##
InstallGlobalFunction(NrOrderedPartitions,function ( arg )
local nr;
if Length(arg) = 1 then
if arg[1] = 0 then
nr := 1;
else
nr := 2^(arg[1]-1);
fi;
elif Length(arg) = 2 then
if arg[1] = 0 then
if arg[2] = 0 then
nr := 1;
else
nr := 0;
fi;
else
nr := Binomial(arg[1]-1,arg[2]-1);
fi;
else
Error("usage: NrOrderedPartitions( <n> [, <k>] )");
fi;
return nr;
end);
#############################################################################
##
#F RestrictedPartitions( <n>, <set> ) . restricted partitions of an integer
##
## 'RestrictedPartitionsA( <n>, <set>, <m>, <part>, <i> )' returns the set
## of all partitions of '<n> + Sum(<part>[[1..<i>-1]])' that contain only
## elements of <set> and that begin with '<part>[[1..<i>-1]]'. To do so it
## finds all elements of <set> that can go at '<part>[<i>]' and calls itself
## recursively for each candidate. <m> is the position of '<part>[<i>-1]'
## in <set>, so the candidates for '<part>[<i>]' are the elements of
## '<set>[[1..<m>]]' that are less than <n>, since we require that
## partitions are nonincreasing.
##
## 'RestrictedPartitionsK( <n>, <set>, <m>, <k>, <part>, <i> )' returns the
## set of all partitions of '<n> + Sum(<part>[[1..<i>-1]])' that contain
## only elements of <set>, that have length '<k>+<i>-1', and that begin with
## '<part>[[1..<i>-1]]'. To do so it finds all elements fo <set> that can
## go at '<part>[<i>]' and calls itself recursively for each candidate. <m>
## is the position of '<part>[<i>-1]' in <set>, so the candidates for
## '<part>[<i>]' are the elements of '<set>[[1..<m>]]' that are less than
## <n>, since we require that partitions are nonincreasing.
##
RestrictedPartitionsA := function ( n, set, m, part, i )
local parts, l;
if n = 0 then
part := ShallowCopy(part);
parts := [ part ];
else
part := ShallowCopy(part);
if n mod set[1] = 0 then
parts := [ part ];
else
parts := [ ];
fi;
for l in [2..m] do
if set[l] <= n then
part[i] := set[l];
Append(parts,RestrictedPartitionsA(n-set[l],set,l,part,i+1));
fi;
od;
if n mod set[1] = 0 then
for l in [i..i+n/set[1]-1] do
part[l] := set[1];
od;
fi;
fi;
return parts;
end;
MakeReadOnlyGlobal( "RestrictedPartitionsA" );
RestrictedPartitionsK := function ( n, set, m, k, part, i )
local parts, l;
if k = 1 then
if n in set then
part := ShallowCopy(part);
part[i] := n;
parts := [ part ];
else
parts := [];
fi;
else
part := ShallowCopy(part);
parts := [ ];
for l in [1..m] do
if set[l]+(k-1)*set[1] <= n and n <= k*set[l] then
part[i] := set[l];
Append(parts,
RestrictedPartitionsK(n-set[l],set,l,k-1,part,i+1) );
fi;
od;
fi;
return parts;
end;
MakeReadOnlyGlobal( "RestrictedPartitionsK" );
InstallGlobalFunction(RestrictedPartitions,function ( arg )
local parts;
if Length(arg) = 2 then
parts := RestrictedPartitionsA(arg[1],arg[2],Length(arg[2]),[],1);
elif Length(arg) = 3 then
if arg[1] = 0 then
if arg[3] = 0 then
parts := [ [ ] ];
else
parts := [ ];
fi;
else
if arg[2] = 0 then
parts := [ ];
else
if not ForAll(arg[2],IsPosInt) then
Error("RestrictedPartitions: Set entries must be positive integers");
fi;
parts := RestrictedPartitionsK(
arg[1], arg[2], Length(arg[2]), arg[3], [], 1 );
fi;
fi;
else
Error("usage: RestrictedPartitions( <n>, <set> [, <k>] )");
fi;
return parts;
end);
#############################################################################
##
#F NrRestrictedPartitions(<n>,<set>) . . . . number of restricted partitions
##
#N 22-Jul-91 martin there should be a better way to do this for given <k>
##
NrRestrictedPartitionsK := function ( n, set, m, k, part, i )
local parts, l;
if k = 1 then
if n in set then
parts := 1;
else
parts := 0;
fi;
else
part := ShallowCopy(part);
parts := 0;
for l in [1..m] do
if set[l]+(k-1)*set[1] <= n and n <= k*set[l] then
part[i] := set[l];
parts := parts +
NrRestrictedPartitionsK(n-set[l],set,l,k-1,part,i+1);
fi;
od;
fi;
return parts;
end;
MakeReadOnlyGlobal( "NrRestrictedPartitionsK" );
InstallGlobalFunction(NrRestrictedPartitions,function ( arg )
local s, n, set, m, p, l;
if Length(arg) = 2 then
n := arg[1];
set := arg[2];
p := [];
for m in [1..n+1] do
if (m-1) mod set[1] = 0 then
p[m] := 1;
else
p[m] := 0;
fi;
od;
for l in set{ [2..Length(set)] } do
for m in [l+1..n+1] do
p[m] := p[m] + p[m-l];
od;
od;
s := p[n+1];
elif Length(arg) = 3 then
if arg[1] = 0 and arg[3] = 0 then
s := 1;
elif arg[1] < arg[3] or arg[3] = 0 then
s := 0;
else
if not ForAll(arg[2],IsPosInt) then
Error("NrRestrictedPartitions: Set entries must be positive integers");
fi;
s := NrRestrictedPartitionsK(
arg[1], arg[2], Length(arg[2]), arg[3], [], 1 );
fi;
else
Error("usage: NrRestrictedPartitions( <n>, <set> [, <k>] )");
fi;
return s;
end);
#############################################################################
##
#F IteratorOfPartitions( <n> )
##
## The partitions of <n> are returned in lexicographic order.
##
## So the partition $\lambda = [ \lambda_1, \lambda_2, \ldots, \lambda_m ]$
## has a successor if and only if $m > 1$.
## If we set $k = \max\{ i; 1 \leq i \leq m-2, \lambda_k > \lambda_{m-1} \}$
## (or $k = 0$ if the set is empty)
## and $l = n - 1 - \sum_{i=1}^{k+1} \lambda_i$
## then the successor of $\lambda$ has the form
## $\mu = [ \lambda_1, \lambda_2, \ldots, \lambda_k, \lambda_{k+1}+1, 1^l ]$
## (where the last term is omitted if $l = 0$).
##
## (Note that $\mu$ is lexicographically larger than $\lambda$,
## clearly $\mu_i = \lambda_i$ for $i \leq k$ is the minimal choice,
## $\mu_{k+1}$ must satisfy $\mu_{k+1} > \lambda_{k+1}$ since
## $\lambda_{k+1} = \lambda_{k+2} = \ldots = \lambda_{m-1} \geq \lambda_m$,
## and for $i > k+1$, $\mu_i = 1$ is the smallest choice.)
##
BindGlobal( "IsDoneIterator_Partitions", iter -> ( iter!.next = false ) );
BindGlobal( "NextIterator_Partitions", function( iter)
local part, m, succ, k;
part:= iter!.next;
m:= Length( part );
if m = 1 then
succ:= false;
else
k:= m-2;
while 0 < k and part[ m-1 ] = part[k] do
k:= k-1;
od;
succ:= part{ [ 1 .. k ] };
k:= k+1;
succ[k]:= part[k] + 1;
Append( succ, 0 * [ 1 .. iter!.n - Sum( succ, 0 ) ] + 1 );
fi;
iter!.next:= succ;
return part;
end );
BindGlobal( "ShallowCopy_Partitions",
iter -> rec( n:= iter!.n, next:= ShallowCopy( iter!.next ) ) );
InstallGlobalFunction( "IteratorOfPartitions", function( n )
if not IsPosInt( n ) then
Error( "<n> must be a positive integer" );
fi;
return IteratorByFunctions( rec(
IsDoneIterator := IsDoneIterator_Partitions,
NextIterator := NextIterator_Partitions,
ShallowCopy := ShallowCopy_Partitions,
n := n,
next := 0 * [ 1 .. n ] + 1 ) );
end );
#############################################################################
##
#F SignPartition( <pi> ) . . . . . . . . . . . . . signum of partition <pi>
##
InstallGlobalFunction(SignPartition,function(pi)
return (-1)^(Sum(pi) - Length(pi));
end);
#############################################################################
##
#F AssociatedPartition( <pi> ) . . . . . . the associated partition of <pi>
##
## 'AssociatedPartition' returns the associated partition of the partition
## <pi> which is obtained by transposing the corresponding Young diagram.
##
InstallGlobalFunction(AssociatedPartition,function(lambda)
local res, k, j;
res := [];
k := Length(lambda);
for j in [1..lambda[1]] do
if j <= lambda[k] then
res[j] := k;
else
k := k-1;
while j > lambda[k] do
k := k-1;
od;
res[j] := k;
fi;
od;
return res;
end);
#############################################################################
##
#F PowerPartition( <pi>, <k> ) . . . . . . . . . . . . power of a partition
##
## 'PowerPartition' returns the partition corresponding to the <k>-th power
## of a permutation with cycle structure <pi>.
##
InstallGlobalFunction(PowerPartition,function(pi, k)
local res, i, d, part;
res:= [];
for part in pi do
d:= GcdInt(part, k);
for i in [1..d] do
Add(res, part/d);
od;
od;
Sort(res);
return Reversed(res);
end);
#############################################################################
##
#F PartitionTuples( <n>, <r> ) . . . . . . . . . <r> partitions with sum <n>
##
## 'PartitionTuples' returns the list of all <r>-tuples of partitions which
## together form a partition of <n>.
##
InstallGlobalFunction(PartitionTuples,function( n, r )
local empty, pm, m, i, s, k, t, t1, res;
empty := rec( tup := List( [1..r], x-> [] ),
pos := List( [1..n-1], x-> 1 ) );
if n = 0 then
return [empty.tup];
fi;
pm := List( [1..n-1], x -> [] );
for m in [ 1 .. QuoInt(n,2) ] do
# the m-cycle in all possible places.
for i in [ 1 .. r ] do
s := rec( tup := List( empty.tup, ShallowCopy ),
pos := ShallowCopy( empty.pos ) );
s.tup[i] := [m];
s.pos[m] := i;
Add( pm[m], s );
od;
# add the m-cycle to everything you know.
for k in [ m+1 .. n-m ] do
for t in pm[k-m] do
for i in [ t.pos[m] .. r ] do
t1 := rec( tup := List( t.tup, ShallowCopy ),
pos := ShallowCopy( t.pos ) );
s := [m];
Append( s, t.tup[i] );
t1.tup[i] := s;
t1.pos[m] := i;
Add( pm[k], t1 );
od;
od;
od;
od;
# collect.
res := [];
for k in [ 1 .. n-1 ] do
for t in pm[n-k] do
for i in [ t.pos[k] .. r ] do
t1 := List( t.tup, ShallowCopy );
s := [k];
Append( s, t.tup[i] );
t1[i] := s;
Add( res, t1 );
od;
od;
od;
for i in [ 1 .. r ] do
s := List( empty.tup, ShallowCopy );
s[i] := [n];
Add( res, s );
od;
return res;
end);
InstallGlobalFunction(NrPartitionTuples, function(n, k)
local res, l, pp, r, a, pr, b;
res := 0;
for l in [0..k] do
pp := Partitions(n, l);
r := Binomial(k, l);
for a in pp do
pr := 1;
for b in a do
pr := pr * NrPartitions(b);
od;
res := res + r * NrArrangements(a, l) * pr;
od;
od;
return res;
end);
#############################################################################
##
#F Lucas(<P>,<Q>,<k>) . . . . . . . . . . . . . . value of a lucas sequence
##
## 'Lucas' uses the following relations to compute the result in $O(log(k))$
## $U_{2k} = U_k V_k, U_{2k+1} = (P U_{2k} + V_{2k}) / 2$ and
## $V_{2k} = V_k^2 - 2 Q^k, V_{2k+1} = ((P^2-4Q) U_{2k} + P V_{2k}) / 2$.
##
InstallGlobalFunction(Lucas,function ( P, Q, k )
local l;
if k = 0 then
l := [ 0, 2, 1 ];
elif k < 0 then
l := Lucas( P, Q, -k );
l := [ -l[1]/l[3], l[2]/l[3], 1/l[3] ];
elif k mod 2 = 0 then
l := Lucas( P, Q, k/2 );
l := [ l[1]*l[2], l[2]^2-2*l[3], l[3]^2 ];
else
l := Lucas( P, Q, k-1 );
l := [ (P*l[1]+l[2])/2, ((P^2-4*Q)*l[1]+P*l[2])/2, Q*l[3] ];
fi;
return l;
end);
##############################################################################
##
#F LucasMod(P,Q,N,k) - return the reduction modulo N of the k'th terms of
## the Lucas Sequences U,V associated to x^2+Px+Q.
##
## Recursive version is a trivial modification of the above function, but
## the running time is dramatically decreased. The running time of the
## the function is dominated by the cost of basic arithmetic operations.
## If reductions mod N are enforced regularly, then these operations are
## constant cost. If not, then they grow quickly as the Lucas sequence
## itself grows exponentially.
##
## See lib/primality.gi for a faster implementation.
##
InstallMethod(LucasMod,
"recursive version, reduce mod N regularly",
[IsInt,IsInt,IsInt,IsInt],
function(P,Q,N,k)
local l;
if k = 0 then
l := [ 0, 2, 1 ];
elif k < 0 then
l := LucasMod( P, Q, N, -k );
if GcdInt( l[3], N ) <> 1 then return fail; fi;
l := [ -l[1]/l[3], l[2]/l[3], 1/l[3] ];
elif k mod 2 = 0 then
l := LucasMod( P, Q, N, k/2 );
l := [ l[1]*l[2], l[2]^2-2*l[3], l[3]^2 ];
else
l := LucasMod( P, Q, N, k-1 );
l := [ (P*l[1]+l[2])/2, ((P^2-4*Q)*l[1]+P*l[2])/2, Q*l[3] ];
fi;
return l mod N;
end);
#############################################################################
##
#F Fibonacci( <n> ) . . . . . . . . . . . . value of the Fibonacci sequence
##
## A recursive Fibonacci needs $O( Fibonacci(n) ) = O(2^n)$ bit operations.
## An iterative version performs $n$ additions, the <i>th involving integers
## with $i$ bits, so we need $\sum_{i=1}^{n}{i} = O(n^2)$ bit operations.
## The binary recursion of 'Lucas' reduces the number of calls to $log2(n)$.
## The number of bit operations now is $O(n)$, i.e., the size of the result.
##
InstallGlobalFunction(Fibonacci,function ( n )
return Lucas( 1, -1, n )[ 1 ];
end);
#############################################################################
##
#F Bernoulli( <n> ) . . . . . . . . . . . . value of the Bernoulli sequence
##
BindGlobal( "Bernoulli2",
[-1/2,1/6,0,-1/30,0,1/42,0,-1/30,0,5/66,0,-691/2730,0,7/6] );
InstallGlobalFunction(Bernoulli,function ( n )
local brn, bin, i, j;
if n < 0 then
Error("Bernoulli: <n> must be nonnegative");
elif n = 0 then
brn := 1;
elif n = 1 then
brn := -1/2;
elif n mod 2 = 1 then
brn := 0;
elif n <= Length(Bernoulli2) then
brn := Bernoulli2[n];
else
for i in [Length(Bernoulli2)+1..n] do
if i mod 2 = 1 then
Bernoulli2[i] := 0;
else
bin := 1;
brn := 1;
for j in [1..i-1] do
bin := (i+2-j)/j * bin;
brn := brn + bin * Bernoulli2[j];
od;
Bernoulli2[i] := - brn / (i+1);
fi;
od;
brn := Bernoulli2[n];
fi;
return brn;
end);
#############################################################################
##
#E combinat.gi . . . . . . . . . . . . . . . . . . . . . . . . . . ends here
##