|
n = 20; (* the value of n is chosen here *)
e = Table[2, {n}]; (*the sequence e*)
Do[
DD = Divisors[k];
e[[k]] = (2^k - Sum[DD[[j]] e[[DD[[j]]]], {j, 1, Length[DD] - 1}])/
k, {k, 2, n}]
PP = IntegerPartitions[n]; npp =
Length[PP]; (*the list of partitions of n*)
(*the maximum length of a cycle in sigma'*)
mlcm = Apply[Max, Table[Apply[LCM, PP[[p]]], {p, npp}]];
(*decompositions of n corresponding to partitions*)
P = Table[0, {i, npp}, {j, n}];
Do[Do[P[[ipp, PP[[ipp, i]]]]++, {i, Length[PP[[ipp]]]}], {ipp, npp}]
EmptyList = Table[0, {j, mlcm}]; (*used to initialize spec(sigma')*)
Vn = 0; Do[(*the main loop through all partitions of n*)
PPP = PP[[p]]; np = Length[PPP]; (*current partition*)
Spec = EmptyList; (*initialization of spec(sigma')*)
divsets = {};
nd = 1;
Do[(*k is the index of the current Partition element*)
DD = Divisors[PPP[[k]]];
AppendTo[divsets, DD];
nd *= Length[DD], {k, 1, np}];
(*divsets is the list of the sets of divisors of cycle lengths in \
sigma*)
Descartes = Tuples[divsets]; (* nd is the length of Descartes *)
Do[ (*loop through Descartes product *)
product = Descartes[[id]];
npr = Length[product];
lcm = 1; prx = 1; pry = 1;
(* Theorem 2 *)
Do[
lcm = LCM[lcm, product[[ipr]]];
prx *= product[[ipr]];
pry *= e[[product[[ipr]]]], {ipr, npr}];
Spec[[lcm]] += prx*pry/lcm, {id, nd}];
(* Theorem 1 *)
numerator = Product[i^Spec[[i]] Spec[[i]]!, {i, Length[Spec]}];
denominatorr = Product[i^P[[p, i]] P[[p, i]]!, {i, n}];
sum = numerator/denominatorr^2;
Vn += sum, {p, npp}]
| 