(* 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)];