upto(n) = {
	my(maxx = 0, fourthrootn = sqrtnint(n, 4), res = List(), c, maxy = sqrtint(n\3)); 
	for(y = 1, maxy,
		maxx = y + fourthrootn; 
		for(x = y + 1, maxx,
			c = (x-y)^2*(x^2 + x*y + y^2);
			
			if(c <= n,
				if(is(c),
					
					listput(res, c)
				);
			,
				break
			)
		)
	);
	Set(res)
}

is(n) = {
	if(valuation(n, 3) == 1, return(0));
	my(f = factor(n), cf = f, q, c, dc);
	cf[,2]>>=1; 
	c = factorback(cf);
	dc = divisors(c); 
	for(i = 1, #dc, 
		dc2 = dc[i]^2;
		dk = n/dc2;
		if(dk > dc2 && (dk - dc2)%3 == 0,
			D = dc2 + 4*(dk - dc2)/3;
			if(issquare(D, &sD) && denominator((-dc[i] + sD)/2) == 1,
				q++
			)
		)	
	);
	q >= 2
	
}