让我们从一个判断一个数是否为素数的函数开始；我们将使用 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 的素数计数来初始化缓存（这是一个糟糕的双关语，抱歉）：
primeCount = [0,0,1]
函数pCount(n)计算 Collatz 序列中n的素数。一旦序列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))
多么有趣的休闲数学题啊！
编辑：由于人们问起米勒-拉宾素性测试仪，让我展示这个基于 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
Sayingfilter(isPrime,range(1000000))在几秒钟内识别出 78498 个少于一百万的素数。您可能希望将此素性测试仪的时序与 Miller-Rabin 测试仪进行比较，并根据运行时效率确定交叉点的位置；我没有做任何正式的计时，但它似乎与 Miller-Rabin 测试仪几乎同时解决了 65-collatz-prime 问题。或者，您可能希望将试验划分与 2、3、5 轮组合到某个限制，例如 1000 或 10000 或 25000，然后对幸存者运行 Miller-Rabin；这可以快速消除大多数复合材料，因此它可以非常快速地运行。