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);