with(numtheory);
with(group):
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_cycleind_edg :=
proc(n)
        option remember;
        local s, t, res, cycs, l1, l2, flat, u, v;

        if n=0 then return 1; fi;

        s := 0:
        for t in pet_cycleind_symm(n) do
            flat := pet_flatten_term(t);

            cycs := flat[2]; res := 1;

            for u to nops(cycs) do
                for v from u+1 to nops(cycs) do
                    l1 := op(1, cycs[u]); l2 := op(1, cycs[v]);
                    res := res * a[lcm(l1, l2)]^(l1*l2/lcm(l1, l2));
                od;
            od;

            for u to nops(cycs) do
                l1 := op(1, cycs[u]);
                if type(l1, odd) then
                   # a[l1]^(1/2*l1*(l1-1)/l1);
                   res := res*a[l1]^(1/2*(l1-1));
                else
                   # a[l1/2]^(l1/2/(l1/2))*a[l1]^(1/2*l1*(l1-2)/l1)
                   res := res*a[l1/2]*a[l1]^(1/2*(l1-2));
                fi;
            od;

            s := s + flat[1]*res;
        od;

        s;
end;



pet_varinto_power :=
proc(N, ZA)
           local beta, res, indvars, v, t, d, m, pot, ZB;

           ZB := pet_cycleind_symm(N);
           if N=1 then ZB := [ZB]; fi;

           indvars := indets(ZA);

           res := 0;
           for beta in ZB do
               t := [];
               for v in indvars do
                   pot := op(1, v);

                   m := 0;
                   for d in divisors(pot) do
                       m := m +
                       d*degree(beta, op(0,v)[d]);
                   od;

                   t := [op(t), v=1+m*z^pot];
               od;

               res := res+lcoeff(beta)*subs(t, ZA);
           od;

           res;
end;

vgf :=
proc(N, q)
         option remember;
         local gf;

         gf := pet_varinto_power(N, pet_cycleind_edg(q));
         expand(gf);
end;

v := (N,q) -> coeff(vgf(N, q), z, q*(q-1)/2);