possiblePrimes(lim,p)={
	my(v=List(),t);
	for(e=0,log(lim)\log(p),
		t=p^e;
		while(t<=lim,
			if(ispseudoprime(t+1) && valuation(t,p)>1, listput(v,t+1));
			t*=2*p
		)
	);
	vecsort(Vec(v))
};
addhelp(possiblePrimes, "possiblePrimes(lim,p): List primes q <= lim+1 for which phi(q) | p^e for some e > 1.");


A053576(n)={
	if(n>31,
		if(n >= 8589934592 && valuation(n>>5, 2)>27,
			warning("Result is conjectural on the nonexistence of Fermat primes >= F(33).")
		);
		return(2<<n)
	);
	n=binary(n);
	prod(i=1,#n,(2^2^(i-1)+1)^n[#n+1-i])
};
addhelp(A053576, "A053576(n): Returns the least number m such that phi(m) = 2^n.");


aPrime(p)={
	my(v=[],u,f=q->q*A053576(valuation(q-1,p)-valuation(q-1,2)),lim=p);
	while(#v==0,
		v=possiblePrimes(lim*=lim,p)
	);
	u=apply(f,v);
	if(vecmin(u)>lim,
		v=possiblePrimes(vecmin(u),p);
		for(i=#u+1,#v,u=concat(u,f(v[#u+1])))
	);
	valuation(eulerphi(vecmin(u)),2)
};
addhelp(aPrime, "aPrime(p): Computes A227533(p) for p an odd prime other than a Fermat prime.");


A227533(n)={
	my(k=2);
	if(isprime(n),
		if(n>257 && n!=65537,
			if(n>>valuation(n-1,2)==1,
				warning("Result is conjectural on the nonexistence of Fermat primes >= F(33).")
			);
			return(aPrime(n))
		)
	);
	while(!istotient((2*n)^k), k++);
	k
};
addhelp(A227533, "A227533(n): Returns the smallest e > 1 such that (2n)^e is a totient, if such e exists; otherwise loop forever. Includes an efficient procedure for handling prime n, which are otherwise slow.");