<p>让我们从一个确定一个数是否为素数的函数开始;我们将使用Miller-Rabin算法,它比我们要处理的数字大小的试除法快:</p>
<pre><code>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
</code></pre>
<p>由于我们想找到满足条件的<em>第一个</em>数,我们的策略是从2开始,逐步向上,并在进行时存储结果。我们用Collatz序列中从0到2的素数来初始化缓存(这是一个糟糕的双关语,抱歉):</p>
^{pr2}$
<p>函数<code>pCount(n)</code>计算Collatz序列中的素数。一旦序列的当前值<em>k</em>下降到<em>n</em>以下,我们就在缓存中查找结果。在此之前,我们测试Collatz序列中的每个奇数是否素数,并在适当的情况下增加素数<em>p</em>。当我们有了<em>n</em>的素数时,我们将其添加到缓存并返回它。在</p>
<pre><code>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
</code></pre>
<p>现在只需计算每个<em>n</em>的素数,从3开始,当素数超过65时停止:</p>
<pre><code>n = 3
t = pCount(n)
while t < 65:
n = n + 1
t = pCount(n)
</code></pre>
<p>用不了多久,在我的电脑上不到一分钟。结果如下:</p>
<pre><code>print n
</code></pre>
<p>结果中有67个素数。如果您想看到它们,这里有一个简单的函数,它打印给定<em>n</em>的Collatz序列:</p>
<pre><code>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
</code></pre>
<p>下面是素数的列表:</p>
<pre><code>filter(isPrime,collatz(n))
</code></pre>
<p>多么有趣的休闲数学问题!在</p>
<p><strong>编辑:</strong>既然有人问过Miller-Rabin素性测试仪,让我展示一下这个基于2,3,5轮的简单素性测试仪;它用2,3,5和不是2,3,5的倍数的数字进行试除法,这包括一些复合物,所以用素数进行试除法效率有点低,但是不需要预先计算和存储素数,所以使用起来更容易。在</p>
<pre><code>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
</code></pre>
<p>说<code>filter(isPrime,range(1000000))</code>可以在几秒钟内识别出小于100万的78498个素数。您可能需要将这个素性测试仪的计时与Miller-Rabin测试仪进行比较,并根据运行时效率确定交叉点的位置;我没有做任何正式的计时,但它似乎在与Miller-Rabin-tester差不多的时间内解决了65 collatz prime问题。或者你可能想把试飞师和2,3,5轮组合在一起,达到一定的限制,比如1000,10000或25000,然后对幸存者进行米勒拉宾测试;这很快消除了大多数合成物,所以它可以很快运行。在</p>