A354386 - OEIS

%I #11 May 30 2022 10:30:32

%S 3,2,337,2633,14143,6108437,373777931

%N a(n) is the first prime that is the start of a sequence of exactly n primes under the map p -> p + A001414(p-1) + A001414(p+1).

%e a(3) = 337 because 337, 337+A001414(336)+A001414(338) = 383, and 383+A001414(382)+A001414(384) = 593 are prime, but 593+A001414(592)+A001414(594) = 660 is not prime, and 337 is the first prime for which this works.

%p spf:= proc(n) option remember; local t; add(t[1]*t[2], t=ifactors(n)[2]) end proc:

%p f:= n -> spf(n-1)+n+spf(n+1):

%p g:= proc(n) option remember;

%p   if not isprime(n) then return 0 fi;

%p   1 + procname(f(n))

%p end proc:

%p V:= Vector(7): count:= 0:

%p p:= 1:

%p while count < 7 do

%p   p:= nextprime(p);

%p   v:= g(p);

%p   if V[v] = 0 then V[v]:= p; count:= count+1 fi;

%p od:

%p convert(V,list);

%t f[1] = 0; f[n_] := Plus @@ Times @@@ FactorInteger[n]; g[n_] := -1 + Length @ NestWhileList[# + f[# - 1] + f[# + 1] &, n, PrimeQ]; seq[len_, max_] := Module[{s = Table[0, {len}], c = 0, p = 1, i}, While[p < max && c < len, p = NextPrime[p]; i = g[p]; If[i <= len && s[[i]] == 0, c++; s[[i]] = p]]; s]; seq[6, 10^7] (* _Amiram Eldar_, May 29 2022 *)

%o (Python)

%o from sympy import factorint, isprime, nextprime

%o def A001414(n): return sum(p*e for p, e in factorint(n).items())

%o def f(p): return p + A001414(p-1) + A001414(p+1)

%o def a(n):

%o     p, count = 1, 0

%o     while count != n:

%o         p = nextprime(p)

%o         fp, count = f(p), 1

%o         while isprime(fp): fp = f(fp); count += 1

%o     return p

%o print([a(n) for n in range(1, 6)]) # _Michael S. Branicky_, May 29 2022

%Y Cf. A001414, A127305.

%K nonn,more

%O 1,1

%A _J. M. Bergot_ and _Robert Israel_, May 24 2022