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