A286658 - OEIS

%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