uphif(f) = prod(i=1, #f~, f[i,1]^f[i,2]-1);
uphin(n) = my(f=factor(n)); uphif(f);



solve_uphi(N, D, limit) = {
	my(g,f,uphi,sol,p,n,pn,uphipn,tmp,ll);
	sol = [] ;
	g = gcd(N, D);
	N/=g;
	D/=g;
	if (D==1,
		if (N==1, sol = [1]);
		sol;
	,
		f = factor(D);
		uphi = prod(i=1, #f~, f[i, 1]^f[i, 2]-1);
		sol = [];	 
		p = f[length(f~),1];
		n = f[length(f~),2];
		pn = p^n;
		uphipn = p^n-1;
		while(pn<=limit,
			tmp = solve_uphi(N*pn, D*uphipn, limit/pn);
			for(i=1,length(tmp),
				if(gcd(pn,tmp[i])==1,
					sol=concat(sol,pn*tmp[i]);
				);
			);
			n++;
			pn*=p;
			uphipn=p^n-1;
		);
		if (uphi == N, sol = concat(sol, [D]));
	);
	sol = vecsort(sol,,8);
	select(x->x<=limit, sol);
}