python - 生成一个列表 a(n) 不是 prime + a(k), k < n 的形式

Question

如何在 Python 中生成这个列表？ a(n) 不是 prime + a(k), k < n 的形式。

这是 oeis 上的列表 http://oeis.org/A025043

它变成 0、1、9、10、25、34、35、49、55、85、91、100、115、121。

我试过大胆的方法，结果并不好。现在我正在寻找一个复杂的解决方案，比如用于素数的埃拉托色尼筛法。 大胆的方法需要迭代每个质数，并且对于质数的每次迭代都需要迭代序列中已经存在的每个数字，这需要很长时间。

这张表是由聪明人生成的:http://oeis.org/A025043/b025043.txt 他们要么使用了大量的计算能力，要么使用了我正在寻找的复杂算法。

To explain what this sequence is, every number that isn't present in it can be represented as a sum of a prime and a number in that sequence. For example 8 = 7(prime) + 1(in sequence), 54 = 53(prime)+1(in sequence), etc.

Answer 1

生成这个序列的明显方法是生成所有素数，然后使用筛子。对于您找到的序列的每个新元素 x，为所需范围内的所有素数 p 删除 x+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 编译提供了一个更快的解决方案。此代码在我的机器上以 7 分 30 秒生成最多 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;
}