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;