with(numtheory);

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;

pet_cycleind_cyclic :=
proc(n)
local s, d;

    s := 0;
    for d in divisors(n) do
        s := s + phi(d)*a[d]^(n/d);
    od;

    s/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;

t :=
proc(n)
option remember;

    if n=1 then return 1 fi;

    1/(n-1)*add(t(n-q)*add(l*t(l), l in divisors(q)),
                q=1..n-1);
end;

U :=
proc(n)
option remember;
local res, m, tgf, dhdgf;

    tgf := add(t(p)*z^p, p=1..n);

    res := 0;

    for m from 3 to n do
        dhdgf :=
        pet_varinto_cind(tgf, pet_cycleind_dihedral(m));
        res := res +
        coeff(expand(dhdgf), z, n);
    od;

    res;
end;