OFFSET
1,1
COMMENTS
Based on a question posed by James Propp. Terms computed by Michael Kleber.
W. Edwin Clark observes (Jun 16 2018) that, based on analysis of the first 10^5 primes, the procedure always ends with {3, 5, 7, 11, 13, 19, 23, 47}, which is sequence A020575. In particular, it appears that the total number of primes obtained is always finite.
First occurrences are in A316226.
REFERENCES
James Propp, Posting to Math Fun Mailing List, Jun 16 2018
LINKS
Hans Havermann, Table of n, a(n) for n = 1..10000
EXAMPLE
a(1)=9: Starting with the first prime, 2, we see that:
2 -> 5 -> 11 -> 23 -> 47 -> 95=5*19,
19 -> 39=3*13,
3 -> 7 -> 15=3*5,
13 -> 27=3*3*3,
which produces 9 different primes, 2 3 5 7 11 13 19 23 47.
MATHEMATICA
propp1[p_] := propp1[p] = #[[1]] & /@ FactorInteger[2*p + 1];
propp[p_Integer] := propp[{p}];
propp[s_List] := propp[s, Union[s, Union @@ propp1 /@ s]];
propp[s_, t_] := If[s == t, s, If[Length[t] > 1000, OVERFLOW[t], propp[t]]];
Table[Length[propp[Prime[n]]], {n, 100}] (* Michael Kleber, Jun 16 2018 *)
g[lst_List] := Union@ Join[lst, First@# & /@ Flatten[FactorInteger[2 lst + 1], 1]]; f[n_] := Length@ NestWhile[g@# &, {Prime@ n}, UnsameQ, All]; Table[ f[n], {n, 100}] (* Robert G. Wilson v, Jun 17 2018 *)
PROG
(Python)
from sympy import prime, primefactors
def a(n):
pn = prime(n)
reach, expand = {pn}, [pn]
while len(expand) > 0:
p = expand.pop()
for q in primefactors(2*p+1):
if q not in reach:
expand.append(q)
reach.add(q)
return len(reach)
print([a(n) for n in range(1, 101)]) # Michael S. Branicky, Jun 29 2022
CROSSREFS
KEYWORD
nonn
AUTHOR
N. J. A. Sloane, Jun 17 2018
STATUS
approved