upto(n) = { my(res = List()); for(e = 2, logint(n, 2), for(d = 2, sqrtnint(n, e), c = d^e + d + e; if(isA254671(c) && isA255265(c) && c <= n, listput(res, c); ); ) ); listsort(res); res } isA255265(n) = { my(d = divisors(n), ind = (#d + 1)\2); while(ind >= 1, D = d[ind]^2 - 4*d[#d + 1 - ind]; if(D < 0, return(0) ); if(issquare(D), return(1) ); ind-- ); 0 } isA254671(n) = { n++; my(qd = numdiv(n)); if(bitand(n, 1), return(qd > 2) , return(qd > 4) ); }