Сито Эратосфена без кратных 2 и 3

Question

Я пытаюсь реализовать алгоритм, описанный в этой статье Faster Sieve of Eratosthenes. Я вполне понимаю эту идею, но не могу понять, как именно она может быть реализована с помощью кода Python.

После некоторой работы я нашел способ конвертировать индексы из сита в сами числа:

number = lambda i: 3 * (i + 2) - 1 - (i + 2) % 2

Но главная проблема - прыжки, которые я должен сделать после получения прайма. Статья объясняет это как:

6np ± p, где p - простое число, а n - некоторое натуральное число.

Есть ли способ описать прыжки с использованием индекса последнего найденного простого числа для такой идеи?

Заранее спасибо.

PS В Objective-C есть реализация. Я довольно новичок в программировании и могу понимать только код на python и js.

Answer 1

Если вы понимаете NumPy, а также Python, посмотрите на эту реализацию primesfrom2to, взятый из этого ответа в Stackru.

def primesfrom2to(n):
    # https://stackru.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Input n>=6, Returns a array of primes, 2 <= p < n """
    sieve = np.ones(n/3 + (n%6==2), dtype=np.bool)
    sieve[0] = False
    for i in xrange(int(n**0.5)/3+1):
        if sieve[i]:
            k=3*i+1|1
            sieve[      ((k*k)/3)      ::2*k] = False
            sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False
    return np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]

В ответе, на который я ссылался, эта процедура была самой быстрой в построении списка простых чисел. Я использую его вариант в своем собственном коде, исследующем простые числа. Детальное объяснение этого заняло бы много места, но оно создает sieve который пропускает кратные 2 и 3. Всего две строки кода (те, которые начинаются sieve[ и окончание = False) выделить кратные числа вновь найденного простого числа. Я думаю, что это то, что вы имеете в виду под "прыжками... после получения прайм" Код хитрый, но изучает работу. Код для Python 2, старый Python.

Вот мой собственный код на Python 3, в котором есть больше комментариев. Вы можете использовать его для сравнения.

def primesieve3rd(n):
    """Return 'a third of a prime sieve' for values less than n that
    are in the form 6k-1 or 6k+1.

    RETURN: A numpy boolean array. A True value means the associated
            number is prime; False means 1 or composite.

    NOTES:  1.  If j is in that form then the index for its primality
                test is j // 3.
            2.  The numbers 2 and 3 are not handled in the sieve.
            3.  This code is based on primesfrom2to in
                <https://stackru.com/questions/2068372/
                fastest-way-to-list-all-primes-below-n-in-python/
                3035188#3035188>
    """
    if n <= 1:  # values returning an empty array
        return np.ones(0, dtype=np.bool_)
    sieve = np.ones(n // 3 + (n % 6 == 2), dtype=np.bool_)  # all True
    sieve[0] = False   # note 1 is not prime
    for i in range(1, (math.ceil(math.sqrt(n)) + 1) // 3): # sometimes large
        if sieve[i]:
            k = 3 * i + 1 | 1  # the associated number for this cell
            # Mark some of the stored multiples of the number as composite
            sieve[k * k                 // 3 :: 2 * k] = False
            # Mark the remaining stored multiples (k times next possible prime)
            sieve[k * (k + 4 - 2*(i&1)) // 3 :: 2 * k] = False
    return sieve

def primesfrom2to(n, sieve=None):
    """Return an array of prime numbers less than n.

    RETURN: A numpty int64 (indices type) array.

    NOTES:  1.  This code is based on primesfrom2to in
                <https://stackru.com/questions/2068372/
                fastest-way-to-list-all-primes-below-n-in-python/
                3035188#3035188>
    """
    if n <= 5:
        return np.array([2, 3], dtype=np.intp)[:max(n - 2, 0)]
    if sieve is None:
        sieve = primesieve3rd(n)
    elif n >= 3 * len(sieve) + 1 | 1:  # the next number to note in the sieve
        raise ValueError('Number out of range of the sieve in '
                         'primesfrom2to')
    return np.r_[2, 3, 3 * np.nonzero(sieve)[0] + 1 | 1]

Спроси меня, есть ли здесь что-то, чего ты не понимаешь.