with(numtheory); P:= proc(q) local a, b, c, n, p;
for n from 2 to q do a:=n*add(op(2, p)/op(1, p), p=ifactors(n)[2]);
b:=ifactors(a)[2]; c:=ifactors(n)[2]; if type(add(c[k][2]/c[k][1], k=1..nops(c))+add(b[k][2]/b[k][1], k=1..nops(b)), integer) then print(n); fi; od; end: P(10^9);