with(numtheory): etr:= proc(p) local b; b:=proc(n) option remember; local d, j; if n=0 then 1 else add(add(d*p(d), d=divisors(j)) *b(n-j), j=1..n)/n fi end end: shr:= proc(p) n->`if`(n=0, 1, p(n-1)) end: b[0]:= etr(n->1): for j from 1 to 7 do b[j]:= etr(shr(b[j-1])) od: a:= shr(b[7]): seq(a(n), n=0..40); # Alois P. Heinz, Sep 08 2008