/*
  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