{ nb = vector(5000, n, numdiv(n^2)); T = vector(#nb, n, vector(nb[n])); occ = vector(#nb); u = 1; i = 0; j = 0; a = vector(10 000); b = 1; for (d = 1, oo, for (k = 1, d, v = lcm(d+1-k, k); i++; if (v <= #T, T[v][occ[v]++] = i; while (occ[u]==nb[u], apply (t -> j++; if (t <= #a, a[t] = j), T[u]); while (a[b], print (b " " a[b]); if (b++ > #a, break (3); ); ); if (u++ > #T, break (3); ); ); ); ); ); } quit