with(combinat);

trees_odd :=
proc(n)
option remember;
local ind, d, a, mset, deg, res;

    if n=1 then return 1 fi;

    res := 0;

    for ind from n^(n-2) to 2*n^(n-2)-1 do
        d := convert(ind, base, n);
        a := [seq(d[q]+1, q=1..n-2)];

        mset := convert(a, `multiset`);
        deg := map(ent -> ent[2]+1, mset);

        if nops(select(q -> type(q, `odd`), deg))
        = nops(deg) then
            res := res + 1;
        fi;
    od;

    res;
end;

X :=
proc(n)
    local q;

    1/2^n*add(binomial(n,q)*(n-2*q)^(n-2), q=0..n);
end;

stirling2even :=
proc(n, k)
    n!*coeftayl(coeftayl(exp(-u+u*1/2*(exp(z)+exp(-z))), u=0, k),
                z=0, n);
end;

stirling2evenA :=
proc(n, k)
    (-1)^(n+k)*k!/2^k*
    add(binomial(n,q)*(-1)^q
        * stirling2(q,k) * stirling2(n-q,k), q=1..n-1);
end;

X2 :=
proc(n)
    add(binomial(n, k)*stirling2evenA(n-2, k)*k!, k=1..n);
end;

X3 :=
proc(n)
    add(binomial(n,k)*(-1)^(n+k)*(k!)^2/2^k*
        add(binomial(n-2,q)*(-1)^q
            * stirling2(q,k) * stirling2(n-2-q,k), q=1..n-3),
        k=1..n);
end;