with(combinat);

pet_cycleind_symm :=
proc(n)
option remember;

    if n=0 then return 1; fi;

    expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;

pet_cycleind_symmNM :=
proc(n, m)
local indA, indB, res, termA, termB, varA, varB,
    lenA, lenB, instA, instB, p, lcmv;
option remember;

    if n=1 then
        return pet_cycleind_symm(m);
    else
        indA := pet_cycleind_symm(n);
    fi;

    if m=1 then
        return pet_cycleind_symm(n);
    else
        indB := pet_cycleind_symm(m);
    fi;

    res := 0;

    for termA in indA do
        for termB in indB do
            p := 1;

            for varA in indets(termA) do
                lenA := op(1, varA);
                instA := degree(termA, varA);

                for varB in indets(termB) do
                    lenB := op(1, varB);
                    instB := degree(termB, varB);

                    lcmv := lcm(lenA, lenB);
                    p :=
                    p*a[lcmv]^(instA*instB*lenA*lenB/lcmv);
                od;
            od;

            res := res + lcoeff(termA)*lcoeff(termB)*p;
        od;
    od;

    res;
end;



mat_count :=
proc(n, m, q)

    subs([seq(a[p]=q, p=1..n*m)],
         pet_cycleind_symmNM(n, m));
end;