with(combinat); with(numtheory);

pet_cycleind_symm :=
proc(n)
option remember;
local l;

    if n=0 then return 1; fi;

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


pet_cycleind_cyclic :=
proc(n)
local d;

    add(phi(d)*a[d]^(n/d), d in divisors(n))/n;
end;

pet_cycleind_dihedral :=
proc(n)
local s;

    s := 1/2*pet_cycleind_cyclic(n);

    if(type(n, odd)) then
        s := s + 1/2*a[1]*a[2]^((n-1)/2);
    else
        s := s + 1/4*(a[1]^2*a[2]^((n-2)/2) + a[2]^(n/2));
    fi;

    s;
end;


pet_prod2vrep :=
proc(varp)
    local v, d, q, res, len, cyc, p;

    q := 1; res := [];

    for v in indets(varp) do
        d := degree(varp, v);
        len := op(1, v);

        for cyc to d do
            res :=
            [op(res), mul(A[p], p=q..q+len-1)];
            q := q+len;
        od;
    od;

    res;
end;

BRACL_CONTRIB :=
proc(gf, n, k)
local q, contr;

    contr := gf;
    for q to n do
        contr := diff(contr, A[q]$k);
        contr := 1/k!*subs(A[q]=0, contr);
    od;

    contr;
end;

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

    if n = 0 or n = 1 then return 1 fi;

    idx_slots := pet_cycleind_dihedral(k*n);
    idx_colors := pet_cycleind_symm(n);

    res := 0;

    for term_b in idx_colors do
        rep := pet_prod2vrep(term_b);

        for term_a in idx_slots 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 cyc_b in rep do
                    len_b := nops(cyc_b);

                    if len_a mod len_b = 0 then
                        if len_a/len_b <= k then
                            q := q + len_b*cyc_b^(len_a/len_b);
                        fi;
                    fi;
                od;

                p := p*q^inst_a;
            od;

            res := res +
            lcoeff(term_a)*lcoeff(term_b)
            *BRACL_CONTRIB(p, n, k);
        od;
    od;

    res;
end;

BFILE :=
proc(mx, k)
local fd, fname, n;

    fname := sprintf("bracl%d-bfile.txt", k);
    fd := fopen(fname, WRITE);

    for n from 0 to mx do
        fprintf(fd, "%d %d\n", n, BRACL(n, k));
    od;
    
    fclose(fd);
end;