with(numtheory);
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_flatten_termA :=
proc(varp)
local terml, d, cf, v;

    terml := [];

    cf := varp;
    for v in indets(varp) do
        d := degree(varp, v);
        terml := [op(terml), [op(1,v), d]];
        cf := cf/v^d;
    od;

    [cf, terml];
end;


EX :=
proc(M,N,T)
option remember;
local res, len_a, inst_a,
    idx_cols, b, flat_b, cycs_b, lcm_b, gf,
    d, f;

    if M*N < T then return 0 fi;
    if T = 1 then return 1 fi;

    idx_cols := pet_cycleind_symm(T);

    res := 0;

    for b in idx_cols do
        flat_b := pet_flatten_termA(b);
        cycs_b := flat_b[2];

        lcm_b := lcm(seq(op(1, v), v in cycs_b));

        for d in divisors(M) do
            for f in divisors(N) do
                len_a := lcm(d, f);
                inst_a := gcd(d, f)*M/d*N/f;

                if len_a mod lcm_b = 0 and
                inst_a >= degree(b)  then
                    gf :=
                    mul((exp(z*op(1, cycs_b[q]))-1)
                        ^op(2, cycs_b[q]),
                        q=1..nops(cycs_b));

                    res := res +
                    1/M/N*phi(d)*phi(f)*flat_b[1] *
                    inst_a! *
                    coeftayl(gf, z=0, inst_a);
                fi;
            od;
        od;
    od;

    res;
end;

ENUM :=
proc(M,N,T)
option remember;
local ind, d, all, orbit, orbitA, orbits, rotM, rotN,
    pos, row, col, perm, conf, pconf;

    if T = 1 then return 1 fi;

    all := M*N;
    orbits := table();

    for ind from T^all to 2*T^all-1 do
        d := convert(ind, base, T)[1..all];

        if nops(convert(d, `set`)) < T then
            next;
        fi;

        orbit := [];

        for rotM from 0 to M-1 do
            for rotN from 0 to N-1 do
                perm := [];
                for pos from 0 to all-1 do
                    col := pos mod N;
                    row := (pos-col)/N;

                    perm :=
                    [op(perm),
                     (row+rotM mod M)*N+
                     (col+rotN mod N)];
                od;

                orbit :=
                [op(orbit),
                 [seq(d[perm[q]+1], q=1..all)]];
            od;
        od;

        orbitA := Array(1..M*N*T!); pos := 1;

        perm := firstperm(T);
        while type(perm, `list`) do
            for conf in orbit do
                pconf :=
                subs([seq(q-1=perm[q]-1, q=1..T)], conf);

                orbitA[pos] := pconf;
                pos := pos + 1;
            od;

            perm := nextperm(perm);
        od;

        orbits[sort(orbitA)[1]] := 1;
    od;

    numelems(orbits);
end;


A := (mx, T) ->
seq(seq(EX(M,N,T), N=1..M), M=1..mx);