with(combinat);

pet_cycleind_symm_maxcl :=
proc(n, m)
option remember;

    if n=0 then return 1; fi;

    expand(1/n*add(a[l]*pet_cycleind_symm_maxcl(n-l, m),
                   l=1..min(m,n)));
end;


pet_cycleind_symm :=
proc(n)
local l;
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_varinto_cind :=
proc(poly, ind)
local subs1, subs2, polyvars, indvars, v, pot, res;

    res := ind;

    polyvars := indets(poly);
    indvars := indets(ind);

    for v in indvars do
        pot := op(1, v);

        subs1 :=
        [seq(polyvars[k]=polyvars[k]^pot,
             k=1..nops(polyvars))];

        subs2 := [v=subs(subs1, poly)];

        res := subs(subs2, res);
    od;

    res;
end;

T1 :=
proc(n, k, m)
option remember;
local rep, q, p, gf;

    rep := -1 + mul(add(A[q]^p, p=0..m), q=1..n);
    gf := pet_varinto_cind(rep, pet_cycleind_symm_maxcl(k, m));
    gf := expand(gf);

    for q to n do
        gf := coeff(gf, A[q], m);
    od;

    gf;
end;

# sanity check for small arguments of the parameters
ENUM :=
proc(n, k, m)
option remember;
local mset, allmsets, idx, digits, dix, src, sidx;

    if k=1 then return 1 fi;

    src := [seq(V[q]$m, q=1..n)];
    allmsets := table();

    for idx from k^(m*n) to 2*k^(m*n)-1 do
        digits := convert(idx, `base`, k)[1..m*n];
        if nops(convert(digits, `set`)) = k then
            mset := table([seq(q=1, q=1..k)]);
            for sidx to m*n do
                dix := digits[sidx] + 1;
                mset[dix] := mset[dix] * src[sidx];
            od;

            allmsets[sort([entries(mset, `nolist`)])] := 1;
        fi;
    od;

    numelems(allmsets);
end;


# best answer to problem
pet_cycleind_symrow :=
(n, m) -> expand(pet_cycleind_symm(m)^n);

T2aux :=
proc(n, k, m)
option remember;
local idx_slots, idx_colors, res, term_a, term_b,
    v_a, v_b, inst_a, inst_b, len_a, len_b, p, q;

    if k = 1 then return 1 fi;

    idx_slots := pet_cycleind_symrow(n,m);
    idx_colors := pet_cycleind_symm(k);

    if m=1 then idx_slots := [idx_slots]; fi;
    
    res := 0;

    for term_a in idx_slots do
        for term_b in idx_colors do
            p := 1;

            for v_a in indets(term_a) do
                len_a := op(1, v_a);
                inst_a := degree(term_a, v_a);

                q := 0;

                for v_b in indets(term_b) do
                    len_b := op(1, v_b);
                    inst_b := degree(term_b, v_b);

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

                p := p*q^inst_a;
            od;

            res := res +
            lcoeff(term_a)*lcoeff(term_b)*p;
        od;
    od;

    res;
end;

T2 :=
proc(n,k,m)
    if k=1 then return 1 fi;
    T2aux(n,k,m)-T2aux(n,k-1,m);
end;


BFILE :=
proc(MXn,m)
    local n, k, idx, fname, fd;

    fname :=
    sprintf("multiset-gen%d-mx%d.bfile", m, MXn);
    
    idx := 1;
    fd := fopen(fname, `WRITE`);
    for n to MXn do
        for k to m*n do
            fprintf(fd, "%d %d\n", idx, T2(n,k,m));
            idx := idx+1;
        od;
    od;
    fclose(fd);
end;