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(A[p]*A[r], r=p+1..2*q-1), p=0..2*q-1)
    - add(A[2*p]*A[2*p+1], 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 :=
proc(q)
option remember;
local res, sl, recurse, edgecycs, 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 to q-1 do
        sl :=
        [seq(A[p]=A[(p+2*rotind) mod 2*q], p=0..2*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;
end;