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;

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;

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

bs_tiles_dh :=
proc(q)
local res, sl, recurse, cind, vsets,
    edges, edgeperm, rotind, seen, factored;

    edges := edges_all_src(q);

    res := tiles_no_symm(q);

    recurse :=
    proc(cycs, sofar, pos)
    local data, nxt;

        if nops(sofar) = 2*q then
            return 1
        fi;

        if pos > nops(cycs) then
            return 0
        fi;

        data := recurse(cycs, sofar, pos+1);

        if nops(cycs[pos] intersect sofar) = 0 then
            data := data +
            recurse(cycs, cycs[pos] union sofar, pos+1);
        fi;

        return data;
    end;

    cind := [];

    seen := table();
    for rotind from 0 to q-1 do
        if rotind >= 1 then
            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)];
        fi;

        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;

    for factored in cind do
        if not(type(seen[factored[2]], `integer`)) then
            vsets :=
            map(indets,
                select(cf -> { seq(degree(cf, v), v in indets(cf)) }
                       = { 1 },
                       map(cyc -> mul(el, el in cyc), factored[1])));

            seen[factored[2]] := recurse(vsets, {}, 1);
        fi;

        res := res + seen[factored[2]];
    od;

    res/q/2;
end;