Z := proc(n) Z(n) := 1/n*expand(add(a[k]*Z(n-k),k=1..n)) end proc;
Z(0) := 1;
magmas := proc(n)
    #option remember;
    uses NumberTheory;
    if n=0 or n=1 then 1 else
    q := 0;
    for m in Z(n) do
        r := 1;
        S := indets(m);
        for u in indets(m) do
            i := op(1,u);
            j := degree(m,u);
            D := 0;
            for d in Divisors(i) do D := D + d*degree(m,a[d]) od;
            r := r*D^(i*j^2);
            S := S minus {u};
            for v in S do
                k := op(1,v);
                l := degree(m,v);
                D := 0;
                for d in Divisors(ilcm(i,k)) do D := D + d*degree(m,a[d]) od;
                r := r*D^(igcd(i,k)*j*l*2);
            od;
        od;
        q := q + coeffs(m)*r;
    od;
    fi;
end proc;