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;

G :=
proc(n,k)
option remember;
local all, term, termvars, res, l1, l2, inst1, u, v,
    uidx, vidx;

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

    all := 0:
    for term in pet_cycleind_symm(n) do
        termvars := indets(term); res := 1;

        # edges on different cycles of different sizes
        for uidx to nops(termvars) do
            u := op(uidx, termvars);
            l1 := op(1, u);

            for vidx from uidx+1 to nops(termvars) do
                v := op(vidx, termvars);
                l2 := op(1, v);

                res := res *
                (k+1)^((l1*l2/lcm(l1, l2))*
                  degree(term, u)*degree(term, v));
            od;
        od;

        # edges on different cycles of the same size
        for u in termvars do
            l1 := op(1, u); inst1 := degree(term, u);
            # a[l1]^(1/2*inst1*(inst1-1)*l1*l1/l1)
            res := res *
            (k+1)^(1/2*inst1*(inst1-1)*l1);
        od;

        # edges on identical cycles of some size
        for u in termvars do
            l1 := op(1, u); inst1 := degree(term, u);
            if type(l1, odd) then
                # a[l1]^(1/2*l1*(l1-1)/l1);
                res := res *
                ((k+1)^(1/2*(l1-1)))^inst1;
            else
                # a[l1/2]^(l1/2/(l1/2))*a[l1]^(1/2*l1*(l1-2)/l1)
                res := res *
                ((k+1)*(k+1)^(1/2*(l1-2)))^inst1;
            fi;
        od;

        all := all + lcoeff(term)*res;
    od;

    all;
end;

C := proc(n, k)
option remember;
local p,q;

    if n < 2 then return 1 fi;

    return G(n,k)
    - 1/n * add(p*C(p,k), p in divisors(n) minus {n})
    - 1/n * add(G(n-q,k)*
                add(p*C(p,k), p in divisors(q)),
                q=1..n-1);
end;