PhiInverse := function(n)
local Solutions;

	Solutions := function(m, primelist)
	local solutions, length, p, e;

	if m = 1 then
		if 2 in primelist then
			return [1, 2];
		else
			return [1];
		fi;
	else
		if Length(primelist) > 1 then
			e := 0;
			p := Last(primelist);
			length := Length(primelist);

			while m mod Phi(p ^ e) = 0 do
				e := e + 1;
			od;

			return Set(Flat(List([0..(e - 1)], i -> Solutions(m / Phi(p ^ i), primelist{[1..(length - 1)]}) * p ^ i)));
		else
			if primelist = [2] and PrimeDivisors(m) = [2] then
				return [2 * m];
			else
				return [];
			fi;
		fi;
	fi;
	end;

return Solutions(n, UnionSet(PrimeDivisors(n), Filtered(DivisorsInt(n), d -> IsPrimeInt(d + 1)) + 1));
end;