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

Question

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

这是 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 您找到的序列，删除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 编译提供了一个更快的解决方案。此代码在我的机器上在 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;
}

Answer 2

Eratosthenes 的筛子看起来与素数筛子非常相似。但是你需要从一个素数列表开始。

对于素数，您有一堆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 的两倍。

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

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