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;