找到给出65个素数的最低collatz序列

Question

我有一个任务，我需要在Python中找到包含65个以上质数的最低Collatz序列。在

例如，19的Collatz序列是：

19, 58, 29, 88, 44, 22, 11, 34, 17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2, 1

7素数序列包含这个。在

我还需要使用回忆录，这样它就不必运行一年来找到它。我找到了Collatz序列记忆的代码，但我不知道当我只需要素数时，如何让它工作。在

这是我找到的Collatz记忆代码：

lookup = {}
def countTerms(n):
    if n not in lookup:
        if n == 1:
            lookup[n] = 1
        elif not n % 2:
            lookup[n] = countTerms(n / 2)[0] + 1
        else:
            lookup[n] = countTerms(n*3 + 1)[0] + 1

    return lookup[n], n

这是我的prime测试仪：

^{pr2}$

Answer 1

让我们从一个确定一个数是否为素数的函数开始；我们将使用Miller-Rabin算法，它比我们要处理的数字大小的试除法快：

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

由于我们想找到满足条件的第一个数，我们的策略是从2开始，逐步向上，并在进行时存储结果。我们用Collatz序列中从0到2的素数来初始化缓存（这是一个糟糕的双关语，抱歉）：

^{pr2}$

函数pCount(n)计算Collatz序列中的素数。一旦序列的当前值k下降到n以下，我们就在缓存中查找结果。在此之前，我们测试Collatz序列中的每个奇数是否素数，并在适当的情况下增加素数p。当我们有了n的素数时，我们将其添加到缓存并返回它。在

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

现在只需计算每个n的素数，从3开始，当素数超过65时停止：

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

用不了多久，在我的电脑上不到一分钟。结果如下：

print n

结果中有67个素数。如果您想看到它们，这里有一个简单的函数，它打印给定n的Collatz序列：

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

下面是素数的列表：

filter(isPrime,collatz(n))

多么有趣的休闲数学问题！在

编辑：既然有人问过Miller-Rabin素性测试仪，让我展示一下这个基于2，3，5轮的简单素性测试仪；它用2，3，5和不是2，3，5的倍数的数字进行试除法，这包括一些复合物，所以用素数进行试除法效率有点低，但是不需要预先计算和存储素数，所以使用起来更容易。在

def isPrime(n): # 2,3,5-wheel
    if n < 2: return False
    wheel = [1,2,2,4,2,4,2,4,6,2,6]
    w, next = 2, 0
    while w * w <= n:
        if n % w == 0: return False
        w = w + wheel[next]
        next = next + 1
        if next > 10: next = 3
    return True

说filter(isPrime,range(1000000))可以在几秒钟内识别出小于100万的78498个素数。您可能需要将这个素性测试仪的计时与Miller-Rabin测试仪进行比较，并根据运行时效率确定交叉点的位置；我没有做任何正式的计时，但它似乎在与Miller-Rabin-tester差不多的时间内解决了65 collatz prime问题。或者你可能想把试飞师和2，3，5轮组合在一起，达到一定的限制，比如1000，10000或25000，然后对幸存者进行米勒拉宾测试；这很快消除了大多数合成物，所以它可以很快运行。在

Answer 2

现有代码缩进错误。我假设这个任务是一个家庭作业任务，所以我不会发布一个完整的解决方案，但我会给你一些有用的片段。在

首先，这里有一个稍微更有效的素性测试程序。它不是测试所有小于a的数字是否都是a的因子，而是测试a的平方根。在

def is_prime(a):
    for i in xrange(2, int(1 + a ** 0.5)):
        if a % i == 0:
            return False
    return True

请注意，此函数为a = 1返回True。没关系，因为您不需要测试1：您可以将它预加载到lookupdict中：

^{pr2}$

您的countTerms函数需要稍作修改，以便在当前n为素数时，它只在lookup计数中添加一个。在Python中，False的数值为0，True的数值为1。这里很方便：

def count_prime_terms(n):
    ''' Count the number of primes terms in a Collatz sequence '''
    if n not in lookup:
        if n % 2:
            next_n = n * 3 + 1
        else:
            next_n = n // 2

        lookup[n] = count_prime_terms(next_n) + is_prime(n)
    return lookup[n]

我已经把函数名改成了Pythonic。在

FWIW，第一个包含65个或更多素数的Collatz序列实际上包含67个素数。它的种子数超过180万个，当检查到该种子的所有序列时，需要进行素性测试的最高数目是151629574372。完成时，lookupdict包含3920492个条目。在

---

为了回应jamesmills关于递归的评论，我写了一个非递归的版本，为了容易地看到迭代版本和递归版本都产生相同的结果，我发布了一个完整的工作程序。我在上面说过我不打算这么做，但我认为现在这样做应该是可以的，因为spørreren已经使用我在原始答案中提供的信息编写了他们的程序。在

我完全同意避免递归是好的，除非它适合于问题域（例如，树遍历）。Python不鼓励递归-它不能优化尾部调用递归，它施加了递归深度限制（尽管如果需要，可以修改该限制）。在

这个Collatz序列素数计算算法自然地递归地表示，但是迭代地做并不难——我们只需要一个列表，在确定序列所有成员的素数时暂时保存序列。这个版本可能需要更多的内存空间。在

当在操作中解决问题时，递归版本达到了343的递归深度。这在默认限制范围内，但仍然不好，如果你想搜索包含更多素数的序列，你会达到这个极限。在

迭代和递归版本的运行速度大致相同（至少在我的机器上是这样）。要解决OP中提到的问题，他们都需要不到2分钟的时间。这比我原来的解决方案快得多，主要是由于素性测试中的优化。在

基本的Collatz序列生成步骤已经需要确定一个数是奇数还是偶数。显然，如果我们已经知道一个数是偶数，那么就没有必要测试它是否是质数。：）我们还可以消除is_prime函数中的偶数因子测试。我们只需将2的结果加载到lookup缓存中，就可以处理2是素数这一事实。在

另一方面，当搜索包含给定数量素数的第一个序列时，我们不需要测试任何从偶数开始的序列。偶数（除了2）不会增加序列的素数，而且由于这样一个序列中的第一个奇数将比我们当前的数字低，所以假设我们从3开始搜索，它的结果将已经在lookup缓存中。如果我们不从3开始搜索，我们只需要确保我们的起始种子足够低，这样我们就不会意外地错过包含所需数量素数的第一个序列。采用这种策略不仅减少了所需的时间，还减少了项目的数量在查找缓存中。在

#!/usr/bin/env python

''' Find the 1st Collatz sequence containing a given number of prime terms

    From http://stackoverflow.com/q/29920691/4014959

    Written by PM 2Ring 2015.04.29

    [Seed == 1805311, prime count == 67]
'''

import sys

def is_prime(a):
    ''' Test if odd `a` >= 3 is prime '''
    for i in xrange(3, int(1 + a ** 0.5), 2):
        if not a % i:
            return 0
    return 1


#Track the highest number generated so far; use a list
# so we don't have to declare it as global...
hi = [2]

#Cache for sequence prime counts. The key is the sequence seed,
# the value is the number of primes in that sequence.
lookup = {1:0, 2:1}

def count_prime_terms_iterative(n):
    ''' Count the number of primes terms in a Collatz sequence 
        Iterative version '''
    seq = []
    while n not in lookup:
        if n > hi[0]:
            hi[0] = n

        if n % 2:
            seq.append((n, is_prime(n)))
            n = n * 3 + 1
        else:
            seq.append((n, 0))
            n = n // 2

    count = lookup[n]
    for n, isprime in reversed(seq):
        count += isprime
        lookup[n] = count

    return count

def count_prime_terms_recursive(n):
    ''' Count the number of primes terms in a Collatz sequence
        Recursive version '''
    if n not in lookup:
        if n > hi[0]:
            hi[0] = n

        if n % 2:
            next_n = n * 3 + 1
            isprime = is_prime(n)
        else:
            next_n = n // 2
            isprime = 0
        lookup[n] = count_prime_terms(next_n) + isprime

    return lookup[n]


def find_seed(numprimes, start):
    ''' Find the seed of the 1st Collatz sequence containing
        `numprimes` primes, starting from odd seed `start` '''
    i = start
    mcount = 0

    print 'seed, prime count, highest term, dict size'
    while mcount < numprimes:
        count = count_prime_terms(i)
        if count > mcount:
            mcount = count
            print i, count, hi[0], len(lookup)
        i += 2


#count_prime_terms = count_prime_terms_recursive
count_prime_terms = count_prime_terms_iterative

def main():
    if len(sys.argv) > 1:
        numprimes = int(sys.argv[1])
    else:
        print 'Usage: %s numprimes [start]' % sys.argv[0]
        exit()

    start = int(sys.argv[2]) if len(sys.argv) > 2 else 3 

    #Round `start` up to the next odd number
    if start % 2 == 0:
        start += 1

    find_seed(numprimes, start)


if __name__ == '__main__':
    main()

当与

$ ./CollatzPrimes.py 65

输出是

seed, prime count, highest term, dict size
3 3 16 8
7 6 52 18
19 7 160 35
27 25 9232 136
97 26 9232 230
171 28 9232 354
231 29 9232 459
487 30 39364 933
763 32 250504 1626
1071 36 250504 2197
4011 37 1276936 8009
6171 43 8153620 12297
10971 44 27114424 21969
17647 48 27114424 35232
47059 50 121012864 94058
99151 51 1570824736 198927
117511 52 2482111348 235686
202471 53 17202377752 405273
260847 55 17202377752 522704
481959 59 24648077896 966011
963919 61 56991483520 1929199
1564063 62 151629574372 3136009
1805311 67 151629574372 3619607