upto(n) = {res = vector(30, i, oo); res[1] = 1; 
for(i = 2, 3, 
	for(j = 1, logint(n, i),
		for(k = 1, max(10000, min(i^j, n \ i^j)),
			m = maxexp(k * i^j); 
			if(m == oo, next);
			if(m > #res, 
				res = concat(res, vector(m - #res, i, oo));  
			); 
			\\print(i" "j" "k" "k * i^j); 
			res[m] = min(res[m], k * i^j);
		)
	) 
); res
}

maxexp(n) = { 
	my(
		fn = factorint(n),
		s = sigma(fn),
		fs = factorint(s)
	);   
	if(denominator(s / factorback(fn[, 1])) != 1,
		return(oo); 
	,
		return(vecmax(vector(#fs~, i, ceil(valuation(n, fs[i, 1])/ fs[i, 2]))))
	)

}