upto(n) = {
	my(res = List(), m = Map(), t = 0);
	for(i = 1, n,
		if(!mapisdefined(m, i),
			if(is(i),
				for(j = 1, sqrtnint(n\i, 2),
					listput(res, i*j^2);
					mapput(m, i*j^2, 1);
				)
			)
		)
	);
	kill(m);
	res
}

is(n) = {
	d = divisors(n);
	d = select(x->x <= sqrtnint(n, 2), d);
	target = n; 
	resforis = 0;
	partialsumd = vector(#d);
	partialsumd[1] = 1;
	for(i = 2, #partialsumd, 
		partialsumd[i] = partialsumd[i-1] + d[i]^2;
	);
	for(i = 1, #d, 
		processforis(d[i]^2, i)	
	);
	resforis
}

processforis(t, ind) = {
	if(resforis == 1,
		return
	);
	if(t == target, 
		resforis = 1;
		return;
	);
	if(target - t > partialsumd[ind],
		return;
	);
	my(ct, cind);
	i = ind - 1;
	for(i = 1, ind-1,
		ct = t + d[i]^2;
		cind = i;
		processforis(ct, cind);
	);
}