Collatz Primes

May 1, 2015

I’ll write the solution in Python because that was the language of the original question on Stack Overflow. Let’s begin with a function that determines if a number is prime; we will use the Miller-Rabin algorithm, which is faster than trial division for the size of numbers that we will be dealing with:

from random import randint

def isPrime(n, k=5): # miller-rabin
    if n < 2: return False
    for p in [2,3,5,7,11,13,17,19,23,29]:
        if n % p == 0: return n == p
    s, d = 0, n-1
    while d % 2 == 0:
        s, d = s+1, d/2
    for i in range(k):
        x = pow(randint(2, n-1), d, n)
        if x == 1 or x == n-1: continue
        for r in range(1, s):
            x = (x * x) % n
            if x == 1: return False
            if x == n-1: break
        else: return False
    return True

Since we want to find the first number that satisfies the criteria, our strategy will be to start from 2 and work our way up, storing results as we go. We prime the cache (that’s a bad pun, sorry) with the prime counts in the Collatz sequences from 0 to 2:

primeCount = [0,0,1]

Function pCount(n) counts the primes in the Collatz sequence for n. As soon as the current value of the sequence k drops below n, we look up the result in the cache. Until then, we test each odd number in the Collatz sequence for primality and increment the prime count p if appropriate. When we have the prime count for n, we add it to the cache and return it.

def pCount(n):
    k, p = n, 0
    while k > 0:
        if k < n:
            t = p + primeCount[k]
            primeCount.append(t)
            return t
        elif k % 2 == 0:
            k = k / 2
        elif isPrime(k):
            p = p + 1
            k = 3*k + 1
        else:
            k = 3*k + 1

Now it’s just a matter of computing the prime counts for each n, in order starting from 3, stopping when the prime count exceeds 65:

n = 3
t = pCount(n)
while t < 65:
    n = n + 1
    t = pCount(n)

That doesn’t take long, less than a minute on my computer. Here’s the result:

>>> print n
1805311

There are 67 primes in the result. If you want to see them, here is a simple function that prints the Collatz sequence for a given n:

def collatz(n):
    cs = []
    while n != 1:
        cs.append(n)
        n = 3*n+1 if n&1 else n//2
    cs.append(1)
    return cs

And here’s the list of primes:

>>> filter(isPrime,collatz(n))
[2707967, 4061951, 6092927, 9139391, 30845447, 46268171, 52051693, 14639539, 1158011, 1737017,
1648811, 2473217, 880361, 660271, 2228417, 352543, 1784753, 47059, 70589, 59561, 100511, 150767,
508841, 381631, 1932011, 2898017, 2173513, 2445203, 2063141, 386839, 580259, 122399, 275399,
348553, 248141, 93053, 117773, 6211, 1747, 2621, 983, 2213, 1579, 1777, 47, 71, 107, 137, 103,
233, 263, 593, 167, 251, 283, 479, 719, 1619, 911, 1367, 577, 433, 61, 23, 53, 5, 2]

You can read more about it at A078350, or run the program at http://ideone.com/Nauozg.

Pages: 1 2

3 Responses to “Collatz Primes”

  1. Rutger said

    In Python..

    Solution performance can be increased by creating a cached collatz(n) function. For example, if we have previously computed collatz(6) already, a consecutive call for collatz(12) would obtain the result from cache.

    def is_prime(number):
        if number == 2:
            return True
        if not (number % 2):
            return False
        i = 3
        while i <= number**0.5:
            if not (number % i):
                return False
            i += 2
        return True
    
    
    def collatz(n):
    	yield n
    	if n == 1:
    		return
    	if n % 2:
    		for c in collatz(3*n+1):
    			yield c
    	else:
    		for c in collatz(n/2):
    			yield c
    
    
    
    for i in range(1,3000000):
    	c = collatz(i)
    	count = 0
    	for v in c:
    		if (is_prime(v)):
    			count +=1
    	if count > 65:
    		print i, count
    

    Result:
    1805311 68
    1993215 71
    2139627 69
    2407081 68
    2707967 68
    2989823 71

  2. Lucas A. Brown said

    In the following Python code, I import an isprime() function from my personal library. This uses BPSW for numbers above 3825123056546413051, Miller-Rabin in various increments for numbers below that, and has a few prefilters including trial division, Fermat’s test, and Euler’s test.

    Brute force:

    from itertools import count
    from labmath import isprime
    def collatz(n):
        yield n
        while n != 1:
            n = 3*n+1 if n%2 else n/2
            yield n
    for n in count(1):
        if sum(map(isprime, list(collatz(n)))) >= 65: break
    print n
    

    On my machine, this finishes in about 15 minutes.

    Caching:

    from itertools import count
    from labmath import isprime
    cache = {1:0}
    for x in count(1):
        n, l = x, []
        while not (n in cache):
            l.append(n)
            n = 3*n+1 if n%2 else n/2
        for m in reversed(l):
            cache[m] = cache[n] + isprime(m)
            n = m
        if cache[x] >= 65: break
    print x
    

    This finishes in 27 seconds.

    (The answer is 1805311, which initiates a sequence of length 402 containing 67 primes.)

  3. Ernie G. said

    I use a sieve to determine if a number, n, is prime where n < 2×10**9. If n is greater than that then I use an isprime method. Since most of the numbers in the Collatz sequence are
    < 2×10**9, the program executes much faster, about 2 sec.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: