(* Mathematica version 4.0.1 *) << DiscreteMath`Combinatorica`; qbinomial[n_Integer, m_Integer] := Expand[Together[Product[(1 - q^(n - j))/(1 - q^(j + 1)), {j, 0, m - 1}]]]; integerpartitions[n_, m_] := Coefficient[qbinomial[n + 1 + m, m], q^n]; rankpartition[(p_)?PartitionQ] := PartitionsP[Tr[p]] - Sum[(integerpartitions[Tr[#1], First[#1] - 1] &)[Drop[p, k]], {k, 0, Length[p] - 1}]; runs[li_List] := Length /@ Split[Sort[li]]; multi[n_, v_, q_, le_] := Take[(Multinomial @@ Append[runs[#1], v - #2] #3) & @@@ ({#, Length[#], If[Length[#] > v || Length[#] < le, 0, 1]} & /@ Partitions[n]), {q, -1}]; multipol[par_?PartitionQ, v_Integer] := multi[Tr[par], v, rankpartition[par], Length[par]]; recur[par_?PartitionQ, v_Integer, h_:{}] := 0 /; v < Length[par]; recur[par_?PartitionQ, v_Integer, h_:{}] := Block[{segs, deco, trabase}, np = Tr[par]; segs = Length /@ Split[TransposePartition[par]]; deco = Table[Max[0, Min[i - j + 1, 1]], {i, 0, #}, {j, #}] & /@ segs; trabase = TransposePartition[par] - 1; Tr[Flatten[ Outer[w[TransposePartition[ DeleteCases[trabase + Flatten[{##}], 0]]] &, Sequence @@ deco, 1]] /. w[q_List] -> temp[q, v - 1, DeleteCases[Sort[Append[h, np - Tr[q] - Tr[h]]], 0]]]]; intermed0[par_List, v_Integer, h_List] := Block[{}, w[Reverse[Sort[Append[#, Tr[h]]]]] & /@ Take[Partitions[ Tr[par]], -(PartitionsP[Tr[par]] + 1 - rankpartition[par])].MapThread[ Times[#1, #2] &, {multipol[par, v], kostka[par]}]]; Clear[kostka]; kostka[{1 ..}] := {1}; kostka[{k_Integer}] := 1 + 0*Range[PartitionsP[k]]; kostka[par_?PartitionQ] := kostka[par] = Block[{npar = Tr[par], v = Tr[par], it}, it = recur[par, v] //. temp[a_, b_, c_] :> recur[a, b, c] /; (Tr[a] == npar && b > 1 || b < Length[a]); vec = ((it /. temp -> intermed0) /. w[p_List] -> w[rankpartition[p]]) + Apply[Plus, w /@ Range[rankpartition[par], PartitionsP[Tr[par]]]]; (((List @@ vec) /. w[_] -> 1) - 1)/(multipol[par, Tr[par]] /. 0 -> 1)];