A316589 - OEIS

%I #21 Sep 07 2018 03:41:31

%S 149,173,307,373,439,443,541,557,563,617,827,863,1297,1303,1373,1453,

%T 1489,1627,1657,1667,1733,1783,1861,1901,2029,2053,2393,2423,2591,

%U 2609,2647,2657,2677,2767,3037,3067,3253,3319,3343,3361,3433,3461,3467,3517,3659

%N Prime numbers p whose number of steps to reach 1 in Collatz (3x+1) problem is a prime number k and, in addition, the least prime number greater than p that also reaches 1 in the same problem in a prime number of steps also does so in k steps.

%e 149 belongs to this sequence as it is prime, it satisfies the Collatz conjecture in 23 (that is prime) steps; no other prime number greater than 149 and less than 163 satisfies the conjecture in a prime number of steps (151 does it in 15 steps; 157 in 36 steps); and the prime number 163 also satisfies it in 23 steps, just as 149 does.

%o (Python)

%o def length_collatz_chain(start):

%o     i=0

%o     while start != 1:

%o         if (start % 2 == 0):

%o             start = start / 2

%o         else:

%o             start = 3 * start + 1

%o         i = i+1

%o     return i

%o def is_prime(num):

%o     if num == 1: return(0)

%o     for k in range(2, num):

%o        if (num % k) == 0:

%o            return(0)

%o     return(1)

%o collatz = []

%o nmax=10000

%o for i in range(nmax):

%o     collatz.append(0)

%o collatz.append(0)

%o for i in range(nmax):

%o     start=i+1

%o     collatz[start]=length_collatz_chain(start)

%o lista_elem=[]

%o elem=[]

%o for i in range(1,nmax):

%o     if is_prime(collatz[i]) and is_prime(i):

%o         elem.append(i)

%o         elem.append(collatz[i])

%o         lista_elem.append(elem)

%o         elem=[]

%o result=""

%o for i in range(len(lista_elem)-1):

%o     if lista_elem[i][1]==lista_elem[i+1][1]:

%o         result=result+str(lista_elem[i][0])+","

%o print(result)

%o (PARI) nbs(n) = my(s); while(n>1, n=if(n%2, 3*n+1, n/2); s++); s; \\ A006577

%o lista(nn) = {vp = primes(nn); vs = select(x->isprime(nbs(x)), vp, 1); vpok = vector(#vs, k, prime(vs[k])); vpoks = vector(#vpok, k, nbs(vpok[k])); for (i=1, #vpoks-1, if (vpoks[i] == vpoks[i+1], print1(vpok[i], ", ")););} \\ _Michel Marcus_, Jul 27 2018

%Y Cf. A006577.

%Y Subsequence of A176112.

%K nonn

%O 1,1

%A _Pierandrea Formusa_, Jul 07 2018