upto(n) = {n = max(120, n); my(t); 
perms = vector(24, i, numtoperm(4, i-1)); 
res = vector(primepi(n)); 
todo = #res - 26; \\ we know 26 primes below 120 aren't of the desired form, so we can stop when we found those. 
termPrimes = primes(primepi(n - (2^4 + 3^3 + 5^2)));
qchecked = 0; 
primeToPrimeIndex = Map(); 

for(i = 1, #termPrimes, 	
	mapput(primeToPrimeIndex, termPrimes[i], i)
);
t = #termPrimes; print(t); 
forprime(p = nextprime(termPrimes[#termPrimes] + 1), n, 
	t++;
	mapput(primeToPrimeIndex, p, t);
);
bounds = concat([[1, 1]], vector(3, i, [2, #termPrimes])); 
forvec(x = bounds, \\ there might be a faster (in terms of finding these primes) "next"-function than the one used by forvec. Maybe by sum of prime indices.  
	for(i = 1, 24, \\ the bounds could probably be set tighter than [1, 24] depending on x and terms not found yet. For example at some point the permutation of exponents
			\\ is mostly if not all the time ending in 1. This might lead to the need of individually verifying a prime isn't there but the need of that seems unlikely
			\\ based on numerical evidence.
		qchecked++; \\ the number we'd like to minimize under the constraint of not missing primes we'd like to find. 
		s = 2^perms[i][1] + sum(j = 2, 4, termPrimes[x[j]] ^ perms[i][j]);
		if(s > n, 
			next(1)
		, 
			
			if(mapisdefined(primeToPrimeIndex, s, &k), \\ see if the number is prime, saving an isprime if it's prime a primepi. 
				if(res[k] == 0,
					todo--;
					res[k] = 1; 
					if(todo%1000 == 0, 
						print(x" "perms[i]" "s" "todo" "qchecked); \\ some progress and an idea of how fast terms are found. 
					); 
					if(todo == 0, return(res)); 
				)
				 
			)
		)
	)
, 
2
); res

}