with(combinat);

pet_cycleind_symm :=
proc(n)
        local p, s;
        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_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;


pet_flat2rep :=
proc(f)
    local p, q, res, cyc, t, len;

    q := 1; res := [];

    for t in f do
        len := op(1, t);
        cyc := [seq(p, p=q+1..q+len-1), q];
        res := [op(res), seq(cyc[p], p=1..nops(cyc))];
        q := q+len;
    od;

    res;
end;

pet_nchoosek_cind :=
proc(n, k)
    option remember;
    local idx_slots, cind, src, aut, q, rep, flat, term;

    cind := 0;

    src := choose(n, k);

    if n=1 then
       idx_slots := [a[1]]
    else
       idx_slots := pet_cycleind_symm(n);
    fi;

    for term in idx_slots do
        flat := pet_flatten_term(term);
        rep := pet_flat2rep(flat[2]);

        aut :=
        map(sel -> sort([seq(rep[sel[q]], q=1..k)]), src);

        cind := cind + flat[1]*pet_autom2cycles(src, aut);
    od;

    cind;
end;


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

    if n=1 then
       idx_cols := [a[1]]
    else
       idx_cols := pet_cycleind_symm(n);
    fi;

    idx_marks := pet_nchoosek_cind(n, k);
    if not type(idx_marks, `+`) then
        idx_marks := [idx_marks];
    fi;

    res := 0;

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

        for b in idx_marks 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;