8 июня 2010 г.

Решето Эратосфена / Sieve of Eratosthenes

Programming Praxis, 19.02.2009

Более двух тысяч лет назад Эратосфен, который вычислил окружность Земли, расстояние до Солнца и наклон земной оси, изобрёл широту и долготу, а также придумал високосный день, создал систематический метод перечисления всех простых числе, который используется и по сей день.
Эратосфен родился в Кирене (на территории современной Ливии) и жил с 276 г. до н. э. по 194 до н. э. Он провёл большую часть своей жизни в Александрии (Египет), где был вторым главой Великой библиотеки перед Аполлонием Родосским.

Решето Эратосфена начинается с создания списка натуральных чисел до желаемого предела; проиллюстрируем метод, вычислив простые числа до тридцати:

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

Теперь возьмём первое число в списке — 2 — и вычеркнем каждое второе число:

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

Затем берём следующее невычеркнутое число в списке — 3 — и вычёркиваем каждое третье число; некоторые из них окажутся уже вычеркнутыми:

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

повторяем последний шаг для следующего невычеркнутого числа — 5:


2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

И так далее, каждый раз вычёркивая все числа, кратные следующему невычеркнутому числу в списке. Все невычеркнутые числа к концу процедуры будут простыми:

2 3 5 7 11 13 17 19 23 29

Этот метод называется решетом, потому что он просеивает весь исходный диапазон чисел в поисках простых чисел.

Решето допускает некоторые оптимизации.

Во-первых, нужно рассматривать только нечётные числа, так как первый проход вычеркнет все чётные за исключением 2, которое можно рассматривать отдельно.

Во-вторых, вычёркивание нужно начинать с квадрата числа, которое рассматривается в данный момент, так как все меньшие составные числа уже были вычеркнуты на предыдущих шагах; например, рассматривая 3, вычёркивание нужно начинать с 9, так как 6 было вычеркнуто, когда рассматривали 2.

В третьих, процедура останавливается на квадратном корне из максимального числа, так как все простые числа, большие него уже будут вычеркнуты к этому моменту. Так, в рассмотренном выше примере не было необходимости проверять числа, большие либо равные 7, так как квадрат 7 больше 30 — наибольшего числа в списке.

Напишите функцию, которая принимает один аргумент n и возвращает список простых чисел, меньших либо равных n используя оптимизированный алгоритм, описанный выше. Примените эту функцию к 15485863 и посчитайте количество полученных простых чисел.


Решение.
Jos Koot написал следующую версию решета Эратосфена:


(define (primes n)
(let* ((max-index (quotient (- n 3) 2))
(v (make-vector (+ 1 max-index) #t)))
(let loop ((i 0) (ps '(2)))
(let ((p (+ i i 3)) (startj (+ (* 2 i i) (* 6 i) 3)))
(cond ((>= (* p p) n)
(let loop ((j i) (ps ps))
(cond ((> j max-index) (reverse ps))
((vector-ref v j)
(loop (+ j 1) (cons (+ j j 3) ps)))
(else (loop (+ j 1) ps)))))
((vector-ref v i)
(let loop ((j startj))
(if (<= j max-index)
(begin (vector-set! v j #f)
(loop (+ j p)))))
(loop (+ 1 i) (cons p ps)))
(else (loop (+ 1 i) ps)))))))


Существует миллион простых числе, меньших, либо равных 15485863:


> (length (primes 15485863))
1000000


Моё решение.

from math import sqrt

def sieve(n):
l = range(3,n+1,2)
i = 0
while 2*i+3 <= sqrt(n):
k = l[i]
l[(k**2-3)/2:(n-1)/2:k] = [0] * ((n-k**2+1)/(2*k)+1)
i += 1
while l[i] == 0:
i += 1
return [2] + filter(lambda x: x != 0, l)

print len(sieve(15485863))

Казалось бы, использование сечений — не самый оптимальный способ присвоить значение сразу многим элементам списка, но при замене страшной строчки «l[(k**2-3)/2:(n-1)/2:k] = [0] * ((n-k**2+1)/(2*k)+1)» на цикл программ стала завершать работу не за 3.172 с, а за 5,056 с.

Всего же для просеивания 15485863 числе программе потребовалось 545 итераций. Из упомянутых трёх секунд только полторы уходило на собственно просеивание. Остальное время программа отфильтровывала вычеркнутые элементы.

Комментариев нет:

Отправить комментарий