Оптимизация алгоритма нахождения простых чисел
2 3 5 7 11 13 17 19 23 29 31… $250.000…
Дело было давно, в университете, когда мы начали изучать язык программирования Pascal и домашним заданием стало создание алгоритма нахождения простых чисел.
Алгоритм был придуман и тутже реализован на изучаемом языке. Программа запрашивала у пользователя число N и искала все простые числа до N включительно. После первого успешного теста сразу же возникло непреодолимое желание ввести N = «много». Программа работала, но не так быстро как хотелось бы. Естественно, дело было в многочисленных проверках (порядка N*N/2), поэтому пришлось избавиться от лишних. В итоге получилось 5 похожих алгоритмов каждый из которых работал быстре предыдущего. Недавно захотелось их вспомнить и реализовать, но на этот раз на Python.
Итак, поехали. Первый алгоритм, ударивший в студенческую голову, продемонстрирован в Листинге 1.
# Листинг 1
# вводим N
n = input("n=")
# создаем пустой список для хранения простых чисел
lst = []
# в k будем хранить количество делителей
k = 0
# пробегаем все числа от 2 до N
for i in xrange(2, n+1):
# пробегаем все числа от 2 до текущего
for j in xrange(2, i):
# ищем количество делителей
if i % j == 0:
k = k + 1
# если делителей нет, добавляем число в список
if k == 0:
lst.append(i)
else:
k = 0
# выводим на экран список
print lst
Очень быстро понимаешь, что в подсчете делителей каждого числа нет никакой надобности и поэтому переменную k можно освободить от своих обязанностей. Действительно, если хотябы один делитель имеется, то число уже не простое. Смотрим Листинг 2.
# Листинг 2
n = input("n=")
lst = []
for i in xrange(2, n+1):
for j in xrange(2, i):
if i % j == 0:
# если делитель найден, число не простое.
break
else:
lst.append(i)
print lst
Конструкция break позволяет нам завершить выполнение внутреннего цикла и перейти к следующей итерации внешнего.
Далее возникает вопрос: «а зачем делить на 4, если на 2 число не делится?». Приходим к выводу, что искать делители нужно только среди простых чисел не превышающих делимое. Наш алгоритм превращается в… см. Листинг 3.
# Листинг 3
n = input("n=")
lst=[]
for i in xrange(2, n+1):
# пробегаем по списку (lst) простых чисел
for j in lst:
if i % j == 0:
break
else:
lst.append(i)
print lst
А потом вспоминаем теорию чисел и понимаем, что переберать надо только числа, не превосходящие корня из искомого. К примеру, если число M имеет делитель pi, то имеется делитель qi, такой, что pi * qi = M. То есть, чтобы найти пару, достаточно найти меньшее. Среди всех пар, предполагаемая пара с максимальным наименьшим — это пара с равными pi и qi, то есть pi * pi = M => pi = sqrt(M). Смотрим Листинг 4.
# Листинг 4
from math import sqrt
n = input("n=")
lst=[]
for i in xrange(2, n+1):
for j in lst:
if j > int((sqrt(i)) + 1):
lst.append(i)
break
if (i % j == 0):
break
else:
lst.append(i)
print lst
Код из Листинга 4 при N=10000 выполняется примерно в 1000 раз быстрее, чем самый первый вариант. Есть еще один «ускоритель», проверять только те числа, которые заканчиваются на 1, 3, 7 или 9 (так как остальные очевидно делятся на 2 или 5). Наблюдаем Листинг 5.
# Листинг 5
from math import sqrt
n = input("n=")
lst=[]
for i in xrange(2, n+1):
if (i > 10):
if (i%2==0) or (i%10==5):
continue
for j in lst:
if j > int((sqrt(i)) + 1):
lst.append(i)
break
if (i % j == 0):
break
else:
lst.append(i)
print lst
В следствии незначительного изменения Листинга 5 получаем небольшую прибавку в скорости:
# Листинг 6
from math import sqrt
n = input("n=")
lst=[2]
for i in xrange(3, n+1, 2):
if (i > 10) and (i%10==5):
continue
for j in lst:
if j > int((sqrt(i)) + 1):
lst.append(i)
break
if (i % j == 0):
break
else:
lst.append(i)
print lst
Итого: Программа из последнего листинга выполняется, примерно, в 1300 раз быстрее первоначального варианта.
Я не ставил перед собой задачи написать программу максимально быстро решающую данную задачу, это скорее демонстрация начинающим программистам того, что правильно составленный алгоритм играет далеко не последнюю роль в оптимизации Ваших программ.
P.S.
Благодаря замечаниям получаем Листинг 7:
# Листинг 7
n = input("n=")
lst=[2]
for i in xrange(3, n+1, 2):
if (i > 10) and (i%10==5):
continue
for j in lst:
if j*j-1 > i:
lst.append(i)
break
if (i % j == 0):
break
else:
lst.append(i)
print lst
при N=10000, поучаем время:
time 1 = 26.24
time 2 = 3.113
time 3 = 0.413
time 4 = 0.096
time 5 = 0.087
time 6 = 0.083
time 7 = 0.053
Решето Эратосфена:
# Листинг 8
n = input("n=")
a = range(n+1)
a[1] = 0
lst = []
i = 2
while i <= n:
if a[i] != 0:
lst.append(a[i])
for j in xrange(i, n+1, i):
a[j] = 0
i += 1
print lst
Результаты при n = 1 000 000:
time 7 = 7.088
time 8 = 1.143
