with(combinat);
with(numtheory);

APBS := n -> add(2^d*mobius(n/d), d in divisors(n));

MXODD_DIV :=
proc(n)
local m;

    m := n;

    while type(m, `even`) do
        m := m/2;
    od;

    m;
end;

APBS_INV := n ->
add(`if`(type(d, `even`), 2^MXODD_DIV(d), 0)*mobius(n/d),
    d in divisors(n));

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

    terml := [];

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

    [cf, terml];
end;

pet_cycleind_bin_p :=
proc(n)
option remember;
local res, prod, term, flat, iter;

    if n=1 then return a[1]^2 fi;

    iter :=
    proc(pos, cycs, sofar)
    local d, len, comb;

        if pos > nops(cycs) then
            len := lcm(seq(q, q in sofar));
            comb := mul(APBS(q), q in sofar);

            prod := prod*a[len]^(comb/len);
            return;
        fi;

        for d in divisors(op(1, cycs[pos])) do
            iter(pos+1, cycs, [op(sofar), d]);
        od;
    end;

    res := 0;

    for term in pet_cycleind_symm(n) do
        flat := pet_flatten_term(term);

        prod := 1;
        iter(1, flat[2], []);
        res := res + flat[1]*prod;
    od;

    res;
end;

bin_p_bool :=
proc(n)
option remember;
local idx, vars;

    idx := pet_cycleind_bin_p(n);
    vars := indets(idx);

    subs([seq(v=2, v in vars)], idx);
end;


pet_cycleind_bin_snp :=
proc(n)
option remember;
local res, prod, term, flat, iter;

    if n=1 then return 1/2*a[1]^2+1/2*a[2] fi;

    iter :=
    proc(pos, cycs, sofar)
    local d, q, len, comb;

        if pos > nops(cycs) then
            len := 1;
            comb := 1;

            for q to nops(cycs) do
                d := op(1, sofar[q]);
                if op(2, sofar[q]) then
                    len :=
                    lcm(len, d/2);
                    comb :=
                    comb *
                    APBS_INV(d);
                else
                    if type(d, `odd`) then
                        len :=
                        lcm(len, 2*d);
                    else
                        len :=
                        lcm(len, d);
                    fi;

                    comb :=
                    comb *
                    (APBS(d)-APBS_INV(d));
                fi;
            od;

            prod := prod*a[len]^(comb/len);
            return;
        fi;

        for d in divisors(op(1, cycs[pos])) do
            if APBS_INV(d) > 0 then
                iter(pos+1, cycs, [op(sofar), [d, true]]);
            fi;
            iter(pos+1, cycs, [op(sofar), [d, false]]);
        od;
    end;

    res := 0;

    for term in pet_cycleind_symm(n) do
        flat := pet_flatten_term(term);

        prod := 1;
        iter(1, flat[2], []);
        res := res + flat[1]*prod;
    od;

    res/2 + pet_cycleind_bin_p(n)/2;
end;

bin_snp_bool :=
proc(n)
option remember;
local idx, vars;

    idx := pet_cycleind_bin_snp(n);
    vars := indets(idx);

    subs([seq(v=2, v in vars)], idx);
end;

bin_snpn_bool :=
proc(n)
option remember;
local idx_slots, idx_cols, res, a, b,
    flat_a, flat_b, cyc_a, cyc_b, len_a, len_b, p, q;

    idx_slots := pet_cycleind_bin_snp(n);
    idx_cols := 1/2*a[1]^2+1/2*a[2];

    res := 0;

    for a in idx_slots do
        flat_a := pet_flatten_term(a);
        for b in idx_cols do
            flat_b := pet_flatten_term(b);

            p := 1;
            for cyc_a in flat_a[2] do
                len_a := op(1, cyc_a);
                q := 0;

                for cyc_b in flat_b[2] do
                    len_b := op(1, cyc_b);

                    if len_a mod len_b = 0 then
                        q := q + len_b;
                    fi;
                od;

                p := p*q;
            od;

            res := res + p*flat_a[1]*flat_b[1];
        od;
    od;

    res;
end;


pet_autom2cycles :=
proc(src, aut)
local numa, numsubs;
local marks, pos, cycs, 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;

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

            while marks[cpos] do
                marks[cpos] := false;
                cpos := numa[cpos];
                clen := clen+1;
            od;

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

        pos := pos+1;
    od;

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

ENUM_cycleind_bin_p :=
proc(n)
option remember;
local cind, perm, src, autom, ind, d;

    src := [];

    for ind from 2^n to 2*2^n-1 do
        d := convert(ind, `base`, 2);
        src := [op(src), d[1..n]];
    od;

    cind := 0;

    perm := firstperm(n);

    while type(perm, `list`) do
        autom :=
        [seq([seq(src[q][perm[p]], p=1..n)],
             q=1..2^n)];

        cind := cind +
        pet_autom2cycles(src, autom);

        perm := nextperm(perm);
    od;

    cind/n!;
end;

ENUM_cycleind_bin_snp :=
proc(n)
option remember;
local cind, perm, src, autom, ind, d;

    src := [];

    for ind from 2^n to 2*2^n-1 do
        d := convert(ind, `base`, 2);
        src := [op(src), d[1..n]];
    od;

    cind := 0;

    perm := firstperm(n);

    while type(perm, `list`) do
        autom :=
        [seq([seq(1-src[q][perm[p]], p=1..n)],
             q=1..2^n)];

        cind := cind +
        pet_autom2cycles(src, autom);

        perm := nextperm(perm);
    od;

    cind/n!/2 + ENUM_cycleind_bin_p(n)/2;
end;