upto(n) = {
	maxr = maxrationOverPhiOfn(n); 
	ulim = n;
	res = List([1]);
	inc = 10^5;
	printat = inc; 
	maxi = sqrtnint(n^2, 3);
	print(maxi);
	
	for(i = 2, maxi,
		if(i > printat,
			print1([i\inc, #res]", ");
			printat += inc;
		);
		process(i)
	);
	Set(res)
}



process(n) = {
	my(f = factor(n));
	if(f[#f~, 1] > maxi,
		return
	);
	f[,2]*=3;
	my(d = divisors(f));
	forstep(i = #d\2, 1, -1,
		if(d[#d + 1 - i] > maxr * d[i],
			return
		);
		if(d[#d + 1 - i] > ulim,
			return
		);
		if(eulerphi(d[#d + 1 - i]) == d[i],
			listput(res, d[#d + 1 - i])
		)
	);
}

maxrationOverPhiOfn(n) = {
	my(res = 1, pp = 1);
	forprime(p = 2, oo, 
		pp*=p;
		if(pp > n, 
			return(ceil(res))
		);	
		res*=(p/(p-1));
	);
}


/* \\ older versions
is(n)=ispower(n*eulerphi(n),3)

upto(n) = {
	inc = 10^6;
	printat = inc;
	res = List();
	forfactored(i = 1, n, 
		if(i[1] > printat, 
			print1(i[1] \ inc", ");
			printat+=inc
		);
		if(ispower(i[1] * eulerphi(i[2]), 3),
			listput(res, i[1]);
			print(); print(res);
		)
	); res
}

upto(n) = {
	res = List();
	forfactored(i = 1, n, 
		if(ispower(i[1] * eulerphi(i[2]), 3),
			listput(res, i[1]);
		)
	); res
}
*/