生成一个列表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。

我试过大胆的方法，结果不太好。现在我正在寻找一个复杂的解决方案，就像伊拉托斯提尼的筛子里的素数。大胆的方法需要对每一个素数进行迭代，而对于每一个素数的迭代，都需要对已经在序列中的每个数进行迭代，这需要很长的时间。

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

为了解释这个序列是什么，它中不存在的每一个数都可以表示为一个素数和一个数列的和。例如，8=7(素数)+1(顺序)，54 =53(素数)+1(顺序)，等等。

Answer 1

生成这个序列的明显方法是生成所有素数，然后使用一个筛子。对于每一个新元素，您所找到的序列的x，为所需范围内的所有素数p剔除x+p。

mathoverflow注释猜测序列的密度趋向于N/sqrt(ln )。如果这是正确的，那么这段代码在时间O(N^2/(ln N)^3/2)中运行。(使用小于N的素数约为N/ln )。

这使得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进行编译，提供了一个更快的解决方案.此代码在我的计算机上生成最多10^7的7m30秒列表：

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

Answer 2

埃拉托斯提尼的筛子看起来和素数的筛子非常相似。但你需要从一个素数的列表开始。

对于素数，您有许多i * prime(k)项，其中prime是我们的序列，i是序列中每个元素的增量。

这里有许多prime(i) + a(k)项，其中a是我们的序列，i是序列中每个元素的增量。

当然，+而不是*在总体效率上有很大的不同。

下面的代码在几秒钟内就能达到100 K，所以如果您想要生成特别大的数字，它会很慢。

你可以等着看是否有人想出了一个更快的方法。

# 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, ...]

作为比较，我写了下面的蛮力方法，一种是迭代素数，检查减法是否给出序列中的元素，另一种是迭代我们的序列，检查我们是否得到一个素数，这两种方法都是100 k的两倍长。

它看起来确实有点像筛子，但它向后看，而不是向前设置值。

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