如何在Python中生成此列表？ a（n）不是素数+ a（k），k <n的形式。
这是oeis http://oeis.org/A025043上的列表
它分为0,1,9,10,25,34,35,49,55,85,91,100,115,121。
我已经尝试过大胆的方式，结果并不顺利。现在我正在寻找一种复杂的解决方案，比如Eratosthenes的Sieve for primes。粗体方式需要迭代每个素数，并且对于素数的每次迭代，迭代序列中已经花费很长时间的每个数。
这个表是由聪明的人生成的：http://oeis.org/A025043/b025043.txt他们要么使用了大量的计算能力，要么使用了复杂的算法，我正在寻找它。
为了解释这个序列是什么，每个不存在的数字可以表示为该序列中的素数和数字之和。例如，8 = 7（素数）+ 1（按顺序），54 = 53（素数）+1（按顺序），等等。
生成此序列的显而易见的方法是生成所有素数，然后使用筛子。对于每个新元素，你找到的序列的x，在所需的范围内击出所有素数x+p的p。
mathoverflow注释猜测序列的密度趋于N / sqrt（ln N）。如果这是正确的，那么该代码在时间O（N ^ 2 /（ln N）^ 3/2）中运行。 （使用小于N的素数的数量约为N / ln N）。
一旦N达到10 ^ 6左右就会变得很慢，但是在我的桌面上运行代码表明即使N = 10 ^ 7，你也会在几个小时内获得完整列表。请注意，算法变得越来越快，所以不要因为最初的缓慢而拖延。
Python 3代码：
import array
N = 10000
def gen_primes(N):
a = array.array('b', [0] * (N+1))
for i in range(2, N+1):
if a[i]: continue
yield i
for j in range(i, N+1, i):
a[j] = 1
def seq(N):
primes = list(gen_primes(N))
a = array.array('b', [0] * (N+1))
for i in range(N+1):
if a[i]: continue
yield i
for p in primes:
if i + p > N:
break
a[i+p] = 1
for i in seq(N):
print(i, end=' ', flush=True)
print()
用C语言重写它，用gcc -O2编译可以提供更快的解决方案。此代码在我的机器上生成7m30s内最多10 ^ 7的列表：
#include <stdio.h>
#include <string.h>
#define N 10000000
unsigned char A[N+1];
int primes[N];
int p_count=0;
int main(int argc, char **argv) {
memset(A, 0, sizeof(A));
for (int i=2; i<=N; i++) {
if(A[i])continue;
primes[p_count++] = i;
for (int j=i; j<=N; j+=i)A[j]=1;
}
memset(A, 0, sizeof(A));
for(int i=0; i<=N; i++) {
if(A[i])continue;
printf("%d ", i);
fflush(stdout);
for(int j=0; j<p_count && i+primes[j]<=N; j++)A[i+primes[j]]=1;
}
return 0;
}
Eratosthenes的Sieve看起来非常类似于素数。但是你需要从一系列素数开始。
对于素数，你有一堆i * prime(k)项，其中prime是我们的序列，i是我们为序列中的每个元素增加的。
在这里你有一堆prime(i) + a(k)术语，其中a是我们的序列，i是我们为序列中的每个元素增加的。
当然，+而不是*在整体效率方面有很大的不同。
下面的代码会在几秒钟内达到100k，所以如果你想生成特别大的数字，它会变慢。
您可以等待，看看是否有人提出了更快的方法。
# taken from https://stackoverflow.com/questions/3939660
def primes_sieve(limit):
a = [True] * limit
a[0] = a[1] = False
for (i, isprime) in enumerate(a):
if isprime:
yield i
for n in range(i*i, limit, i):
a[n] = False
def a_sieve(limit, primes):
a = [True] * limit
for (i, in_seq) in enumerate(a):
if in_seq:
yield i
for n in primes:
if i+n >= limit:
break
a[i+n] = False
primes = list(primes_sieve(100000))
a = list(a_sieve(100000, primes))
# a = [0, 1, 9, 10, 25, 34, 35, 49, 55, 85, 91, 100, 115, 121, 125, 133, 145, ...]
为了比较，我写了下面的暴力方法，一个迭代过质子并检查是否减去给出了我们序列中的一个元素，另一个迭代我们的序列并检查我们是否得到一个素数，两者都需要大约两倍100k。
它看起来有点类似于Sieve，但它向后看，而不是向前设置值。
def a_brute(limit, primes):
a = [True] * limit
for i in range(limit):
for p in primes:
if i-p < 0:
yield i
break
if a[i-p]:
a[i] = False
break
else:
yield i
def a_brute2(limit, primes_sieve):
a = [True] * limit
l = []
for i in range(limit):
for j in l:
if i-j < 0:
l.append(i)
break
if primes_sieve[i-j]:
a[i] = False
break
else:
l.append(i)
return l
结果：
Primes sieve: 0.02 seconds
Sieve: 6.26 seconds
Brute force 1: 14.14 seconds
Brute force 2: 12.34 seconds