%I #19 Jun 16 2022 14:27:15
%S 13,29,53,149,173,269,293,317,389,509,557,653,769,773,797,1109,1229,
%T 1493,1637,1733,1949,1997,2309,2477,2693,2837,2909,2957,3329,3413,
%U 3533,3677,3989,4133,4157,4253,4349,4373,4493,4517,5189,5309,5693,5717,5813,6173
%N Primes of the form p*b^b + 1, where p is a prime and b>1.
%H Robert Israel, <a href="/A286658/b286658.txt">Table of n, a(n) for n = 1..10000</a>
%e a(1) = 3*(2^2)+1 = 13.
%e a(2) = 7*(2^2)+1 = 29.
%e a(3) = 13*(2^2)+1 = 53.
%p N:= 10000: # for all terms <= N
%p Res:= NULL:
%p P:= select(isprime, [2,seq(i,i=3..N/4,2)]):
%p for b from 2 do
%p q:= b^b; if q > N/2 then break fi;
%p for i from 1 to nops(P) do
%p x:= P[i]*q+1;
%p if x > N then break fi;
%p if isprime(x) then Res:= Res, x fi;
%p od od:
%p sort(convert({Res},list)); # _Robert Israel_, Nov 12 2019
%t nmax=10^4; pimax=PrimePi[nmax]; bmax=1;While[(bmax+1)^(bmax+1)<=nmax,bmax++];Select[Union@Flatten@Table[Prime[pi] b^b+1,{b,2,bmax},{pi,pimax}],PrimeQ[#]&&#<=nmax&]
%o (PARI) list(lim)=my(v=List()); lim\=1; for(b=2,oo, my(p=2*b^b+1); if(p>lim, break); if(isprime(p), listput(v,p))); forstep(b=2,oo,2, my(B=b^b); if(3*B+1>lim, break); forprime(q=3,(lim-1)\B, my(p=q*B+1); if(isprime(p), listput(v,p)))); Set(v) \\ _Charles R Greathouse IV_, Jun 16 2022
%Y Cf. A175768, A285015.
%K nonn,easy
%O 1,1
%A _Vincenzo Librandi_, May 12 2017