upto(n) = {
	my(p = 2);
	inc = 10^8;
	printat = inc;
	forprime(q = 3, n,
		if(q > printat,
			print(q \ inc" "resforupto);
			printat += inc;
		);
		if(q - p >= leastunseen,
			c = numberofsquarefreebetweenpandnextprime(p, q);
			if(c > 0,
				if(c > #resforupto,
					resforupto = concat(resforupto, vector(c - #resforupto, i, oo));
				);
				if(p < resforupto[c],
					updateleastunseen();
					resforupto[c] = p
				);
				
			);
		);
		p = q
	);
	resforupto
	
}

updateleastunseen() = {
	my(i);
	for(i = 1, #resforupto,
		if(resforupto[i] == oo,
			leastunseen = i;
			return
		);
	);
	leastunseen = #res + 1
}

numberofsquarefreebetweenpandnextprime(p, q) = {
	if(maxnumberofsquarefreebetween(p + 1, q - 1) < leastunseen,	
		return(0);
	);
	my(res = sum(i = p + 1, q-1, issquarefree(i)));
	res
	
}

maxnumberofsquarefreebetween(n1, n2) = {
	my(res = 0);
	d = divisors(30);
	for(i = 1, #d, 
		res += moebius(d[i]) * max(0, (n2\d[i]^2 - nextmultiple(n1, d[i]^2)\d[i]^2) + 1)
	); res
}

nextmultiple(n, m) = {
	m * (ceil((n)/m))
}