MacMahonsMultisetFormula := proc(L::Or(list,multiset),k::integer)

# Start: 24.10.2010
# Last change: 5.11.2010
# thomas.wieder@t-online.de

# Calculate the number of k-combinations ("selections") from a multiset according to Percy A. MacMahon's formula.
#
# 	   Nmltst
# 	  -----		 -----
# 	  \	    pss	 \
# Result = | 	(-1)	  |	binomial( Nmltst + k - ListOfMultiplicities(i_1) - ListOfMultiplicities(i_2) - ... - ListOfMultiplicities(i_pss) - pss - 1, Nmltst - 1 )
# 	  /		 /
# 	  -----		 -----
# 	  pss=0		 1 <= i_1 < i_2 < ... < i_pss <= Nmltst
#	
# This formula is from Percy Alexander MacMahon.
# Unfortunately, the exact citation is not know to me, please let me know!

# Usage:
#
#> with(combstruct);
#> read "C:/Users/Wieder/Documents/Projekte/multiset/MacMahonsMultisetFormula.mpl";
#> 
#> L := [a, a, a, b, b, b, b, c]: k := 6: MacMahonsMultisetFormula(L, k);
#                                      5
#
# For the example above the five 6-combinations are:
# [a, a, a, b, b, b], [a, a, a, b, b, c], [a, a, b, b, b, b], [a, a, b, b, b, c], [a, b, b, b, b, c].

# Acknowledgement:
# Many thanks to Harun Siljak, University of Sarajevo, Bosnia and Herzegovina, for a helpful discussion.

local AllSelectionsOfIndices, arg1, arg2, ASelectionOfIndices, 
ASelectionOfMultiplicities, iterator, iverbose, ListOfMultiplicities, 
Lmltst, Nmltst, pss, Result, SumForASelectionOfMultiplicities, 
SumOfSelectedMultiplicities; 

iverbose := 0:
# Convert the list into a Maple multiset.
if type(L,list) then
Lmltst := convert(L,multiset);
else
Lmltst := L;
end if;
Nmltst := nops(Lmltst);
if iverbose = 1 then print("Lmltst,Nmltst:",Lmltst,Nmltst); end if;
# Form the list of multiplicities.
ListOfMultiplicities := seq(op(2,Lmltst[i]),i=1..Nmltst);
ListOfMultiplicities := [ListOfMultiplicities];
if iverbose = 3 then print("ListOfMultiplicities:",ListOfMultiplicities); end if;
Result := 0;

# Loop over length of subselection taken from the list of multiplicities.
for pss from 0 to Nmltst do

# The list of multiplicities is a multiset itself.
# We have to form submultisets of the list of multiplicities.
# Maple's commands "choose" and "Combination" do not work on multisets as needed here.
# Thus we have to apply the "Combination" command on the indices of the elements of a multiset, instead on the elements itself.

SumForASelectionOfMultiplicities := 0;

# Loop over the selections with pss parts which are taken from the multiplicites.
iterator := iterstructs(Combination(Nmltst), size=pss);
while not finished(iterator) do 
ASelectionOfIndices := nextstruct(iterator); 
ASelectionOfMultiplicities := [seq(ListOfMultiplicities[ASelectionOfIndices[i]],i=1..pss)];

SumOfSelectedMultiplicities := add(ASelectionOfMultiplicities[i],i=1..pss);
if iverbose = 1 then print("ASelectionOfMultiplicities, SumOfSelectedMultiplicities:",ASelectionOfMultiplicities,SumOfSelectedMultiplicities);  end if;

arg1 := Nmltst + k - SumOfSelectedMultiplicities - pss - 1;
arg2 := Nmltst - 1;
if arg1 > 0 then SumForASelectionOfMultiplicities := SumForASelectionOfMultiplicities + binomial(arg1,arg2); end if;
if iverbose = 1 then print("SumForASelectionOfMultiplicities , arg1, arg2, binomial(arg1,arg2):",SumForASelectionOfMultiplicities , arg1, arg2, binomial(arg1,arg2));  end if;

# End of loop over ASelectionOfMultiplicities.
end do;

Result := Result + ((-1)^(pss)) * SumForASelectionOfMultiplicities;
if iverbose = 1 then print("Result=",Result); end if;
if iverbose = 2 then print("Result, SumForASelectionOfMultiplicities:",Result,SumForASelectionOfMultiplicities); end if;

# End of loop over pss.
end do;

return Result;

end proc;