/* To change the number of terms produced, change these three lines: sf = select(issquarefree, [1..111 000]); u = 1; a = vector(10 000); Replace 10 000 by the number of desired terms; then 111 000 must be greater than the largest term. To get 1000 terms, for example, replace 10 000 by 1000, and 111 000 by the 1000-th prime. */ s = 0 seen(v) = bittest(s, v) see(v) = s = bitor(s, 2^v) { sf = select(issquarefree, [1..111 000]); u = 1; a = vector(10 000); x = 1; for (n=1, #a, v=0; for (i=u, #sf, if (!seen(sf[i]) && gcd(x, sf[i])==1, v=sf[i]; break; ); ); if (v, see(a[n]=v); x*=a[n]; print (n " " a[n]); while (seen(sf[u]), u++; ), break; ); if (n%2==0, x/=a[n/2]; ); ); } quit