with(combinat);

pet_autom2cyclesA :=
proc(src, aut)
local numa, numsubs, marks, pos, cycs,
    data, item, cpos, clen;

    numsubs := [seq(src[k]=k, k=1..nops(src))];
    numa := subs(numsubs, aut);

    marks := Array([seq(true, pos=1..nops(aut))]);

    cycs := []; pos := 1; data := [];

    while pos <= nops(aut) do
        if marks[pos] then
            clen := 0; item := []; cpos := pos;

            while marks[cpos] do
                marks[cpos] := false;
                item := [op(item), aut[cpos]];

                cpos := numa[cpos];
                clen := clen+1;
            od;

            cycs := [op(cycs), clen];
            data := [op(data), item];
        fi;

        pos := pos+1;
    od;

    return [data, mul(a[cycs[k]], k=1..nops(cycs))];
end;


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;


edges_all_src :=
proc(q)
option remember;
local edges;

    edges :=
    add(add(C[p]*C[r], r=p+1..2*q-1), p=0..2*q-1)
    - add(C[2*p]*C[2*p+1], p=0..q-1);

    subs([seq(C[2*p]=A[p], p=0..q-1),
          seq(C[2*p+1]=B[p], p=0..q-1)],
         [seq(p, p in edges)]);
end;

pet_perms_edge_cind :=
proc(q)
option remember;
local cind, perm, sl,
    edges, edgeperm, rotind;

    edges := edges_all_src(q);

    cind := [];

    for rotind from 0 to q-1 do
        sl :=
        [seq(A[p]=A[(p+rotind) mod q], p=0..q-1),
         seq(B[p]=B[(p+rotind) mod q], p=0..q-1)];

        edgeperm := subs(sl, edges);

        cind :=
        [op(cind), pet_autom2cyclesA(edges, edgeperm)];
        
        sl :=
        [seq(A[p]=B[q-1-((p+rotind) mod q)], p=0..q-1),
         seq(B[p]=A[q-1-((p+rotind) mod q)], p=0..q-1)];

        edgeperm := subs(sl, edges);

        cind :=
        [op(cind), pet_autom2cyclesA(edges, edgeperm)];
    od;

    cind;
end;

tiles_dh :=
proc(q)
option remember;
local slot_idx, edge_idx, slot_term, edge_term,
    cycs, deg, gf, var, res;

    slot_idx := `if`(q=1, [a[1]], pet_cycleind_symm(q));
    edge_idx := pet_perms_edge_cind(q);

    res := 0;

    for slot_term in slot_idx do
        for edge_term in edge_idx do
            if type(edge_term[2]/slot_term,
                    `monomial`) then
                gf := 1;

                for var in indets(slot_term) do
                    cycs :=
                    select(c->nops(c)=op(1,var),
                           edge_term[1]);

                    deg := degree(slot_term, var);

                    gf := gf*deg!*op(1, var)^deg
                    *mul(1+z[op(1, var)]*mul(el, el in cyc),
                         cyc in cycs);
                od;

                for var from 0 to q-1 do
                    gf :=
                    subs(A[var] = 0, diff(gf, A[var]));
                    gf :=
                    subs(B[var] = 0, diff(gf, B[var]));
                od;

                for var in indets(slot_term) do
                    gf :=
                    coeftayl(gf, z[op(1, var)]=0,
                             degree(slot_term, var));
                od;

                res := res + lcoeff(slot_term)*gf;
            fi;
        od;
    od;

    res/nops(edge_idx);
end;

tiles_all :=
q -> add(binomial(q,k)*(-1)^k*(2*q-2*k)!/2^(q-k)/(q-k)!,
         k = 0..q);

enum_tiles_all_dh :=
proc(q)
option remember;
local recurse, count, edges, orbits, sl;

    edges := edges_all_src(q);
    orbits := table();

    sl :=
    [seq([seq(A[p]=A[(p+r) mod q], p=0..q-1),
          seq(B[p]=B[(p+r) mod q], p=0..q-1)],
         r=0..q-1)];

    sl :=
    [op(sl),
     seq([seq(A[p]=B[q-1-((p+r) mod q)], p=0..q-1),
          seq(B[p]=A[q-1-((p+r) mod q)], p=0..q-1)],
         r=0..q-1)];     
    
    recurse :=
    proc(sofar, sel, n, pos)
    local orbit, symind;

        if n = q then
            count := count + 1;

            orbit := [];

            for symind to 2*q do
                orbit :=
                [op(orbit),
                 subs(sl[symind], sel)];
            od;

            orbits[sort(orbit)[1]] := 1;

            return;
        fi;

        if pos > nops(edges) then
            return;
        fi;

        recurse(sofar, sel, n, pos+1);

        if nops(indets(sofar*edges[pos])) =
        2*(n+1) then
            recurse(sofar*edges[pos],
                    {op(sel), edges[pos]},
                    n+1, pos+1);
        fi;
    end;

    count := 0;
    recurse(1, {}, 0, 1);

    [count, numelems(orbits)];
end;