big = 1 000 000
s = 0
S = []
unseen = 1
seen(v) = if (v < big, bittest(s, v), setsearch(S, v))
see(v) = if (v < big, s = bitor(s, 2^v), S = setunion(S, [v])); while (seen(unseen), unseen++)

rad(n) = vecprod(factor(n)[,1]~)
ok(r) = {
	forprime (p = 2, oo,
		if (r == 1,
				return (1),
			r % p == 0,
				r /= p,
				return (0);
		);
	);
}
step(ij) = {
	my (s = 1);
	forprime (p = 2, oo,
		if (ij == 1,
				return (s),
			ij % p == 0,
				ij /= p^valuation(ij, p),
				s *= p;
		);
	);
}

other(i,j) = {
	my (s = step(i*j));
	forstep (k = ceil(unseen/s)*s, oo, s,
		if (!seen(k) && gcd(j, k)==1 && ok(rad(i*j*k)),
			return (k);
		);
	);
}

{
	for (n = 1, 10 000,
		if (n <= 2,
				v = unseen,
				v = other(pp, p);
		);

		see(v);

		print (n " " v);

		[pp, p] = [p, v];
	);
}

quit