# Maple programs for A007916 and A278028 (power tower for n),
# A278029, A052409, A089723, A277564, for all n <= N2
# N. J. A. Sloane, Nov 10 2016

# A7916 = list of NPPs (non-prime powers)
# C[n] = height of tower for n; B[n,*] = the tower
# E[n] means n is an E-th power

NN:=30000; N2:=20000;
A7916:=[];
A7916i:=Array(1..2*NN,0);
C:=Array(1..NN,0); E:=Array(1..NN,1);
B:=Array(1..NN,1..4,0);
C[1]:=1;
B[1,1]:=1;
nNPP:=0;
for n from 2 to N2 do
t1:=ifactors(n)[2]; t2:=nops(t1);
g := t1[1][2]; for j from 2 to t2 do g:=gcd(g,t1[j][2]); od:
# do we have a new NPP?
if g=1 then nNPP:=nNPP+1; A7916:=[op(A7916),n]; A7916i[n]:=nNPP; 
C[n]:=1; B[n,1]:=nNPP; 
else 
n2:=simplify(n^(1/g)); E[n]:=g;
   for i from 1 to C[n2] do B[n,i]:=B[n2,i]; od:
   for j from 1 to C[g] do B[n,C[n2]+j]:=B[g,j]; od:
   C[n]:=C[n2]+C[g];
fi;
od:  #od n 
lprint("NPPs");
[seq(A7916[k],k=1..nops(A7916))];  # A007916
lprint("inverse");
[seq(A7916i[k],k=1..100)];         # A278029
lprint("powers");
[seq(E[k],k=1..N2)];               # A052409
lprint("heights");
[seq(C[k],k=1..N2)];               # A089723
lprint("towers");
for n from 1 to 100 do             # A278028 (also A277564)
lprint(n, [seq(B[n,i],i=1..C[n])]); od: