生成列表a(n)不是素数+ a(k),k <n的形式

Question

如何在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(按顺序),等等.

Answer 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;
}