Сформировать список a (n) не имеет вида штрих + a (k), k <n

Question

Как сгенерировать этот список в Python? a (n) не имеет вида простого числа + a (k), k

Вот список на oeis http://oeis.org/A025043

Это идет как 0, 1, 9, 10, 25, 34, 35, 49, 55, 85, 91, 100, 115, 121.

Я попробовал смелый способ, не получилось хорошо. Сейчас я ищу сложное решение, такое как Sieve of Eratosthenes для простых чисел. Жирный способ требует итерации по каждому простому числу, а для каждой итерации простого - итерации по каждому числу уже в последовательности, что занимает очень много времени.

Эта таблица была сгенерирована кем-то умным: 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 дает гораздо более быстрое решение. Этот код генерирует список до 10 ^ 7 в 7m30s на моем компьютере:

#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 000 за несколько секунд, поэтому он будет медленным, если вы хотите генерировать особенно большие числа.

Вы можете подождать и посмотреть, придумает ли кто-нибудь более быстрый метод.

# 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 тыс.

Он выглядит несколько похожим на сито, но смотрит назад, а не устанавливает значения вперед.

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