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_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;

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

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

    for q to n do
        gf := diff(gf, A[q]$m);
        gf := 1/m!*subs(A[q]=0, gf);
    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, q;

    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;

                if src[sidx] in indets(mset[dix]) then
                    break;
                fi;

                mset[dix] := mset[dix] * src[sidx];
            od;

            if sidx = m*n+1 then
                allmsets[sort([entries(mset, `nolist`)])] := 1;
            fi;
        fi;
    od;

    numelems(allmsets);
end;

dgseq :=
proc(ms)
local vars, var, ds, q;
    vars := indets(ms);
    ds := sort([seq(degree(ms, var), var in vars)]);
    mul(A[q]^ds[q], q=1..nops(vars));
end;


T2gen :=
proc(ms, k)
option remember;
local vars, var, rep, term, res, sbms, ell;

    vars := indets(ms);

    if nops(vars) = 0 and k = 0 then
        return 1;
    elif nops(vars) = 0 or k=0 then
        return 0;
    fi;

    if nops(vars) = 1 then
        if degree(ms) = k then
            return 1;
        else
            return 0;
        fi;
    end;

    rep := -1 + expand(mul(1+var, var in vars));
    res := 0;

    for ell to k do
        for term in rep do
            sbms := ms/term^ell;
            if type(sbms, `polynom`) then
                res := res + T2gen(dgseq(sbms), k-ell);
            fi;
        od;
    od;

    res/k;
end;

T2 :=
proc(n,k,m)
    local q;
    T2gen(mul(A[q]^m, q=1..n), k);
end;

OEIS :=
proc(N,m)
local n, k;
    seq(seq(T2(n,k,m), k=m..n*m), n=1..N);
end;

OEISrowsums :=
proc(N,m)
local n, k;
    seq(add(T2(n,k,m), k=m..n*m), n=1..N);
end;


OEIS2 :=
proc(N,m)
local n, k;

    for n to N do
        print(seq(T2(n,k,m), k=m..n*m));
    od;
end;


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

    fname :=
    sprintf("multiset-set%d-mx%d.bfile", m, MXn);

    idx := 1;
    fd := fopen(fname, `WRITE`);
    for n to MXn do
        for k from m to m*n do
            fprintf(fd, "%d %d\n", idx, T2(n,k,m));
            idx := idx+1;
        od;
    od;
    fclose(fd);
end;