seen = 0
unseen = 1

prime powers(n) = my (f=factor(n)); Set(vector(#f~, i, f[i,1]^f[i,2]))

other(p) = {
	seen += 2^p;
	while (bit test(seen, unseen),
		unseen++;
	);
	my (pp=prime powers(p));
	for (v=unseen, oo,
		if (!bit test(seen, v) && set intersect(prime powers(v), pp),
			return (v);
		);
	);
}

for (n=1, 10 000, v = if (n<=2, n, other(v)); print (n " " v))

quit