Алгоритм нахождения простых чисел

Оптимизация алгоритма нахождения простых чисел


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

Similar posts

Ads
AdBlock has stolen the banner, but banners are not teeth — they will be back

More

Comments 64

    +16
    Сделайте листинг 7, в листинге 6 вынесите
    int((sqrt(i)) + 1)
    

    За пределы внутреннего цикла, т.к. вы впустую вызываете очень тяжёлую функцию корня. Другой вариант — это не вызывать тяжёлую функцию корня, а хранить ещё один массив из квадратов простых чисел, сравнивать размер с которым. Как-то так:
    if (j * j - 1) > i:
    

    Должно стать значительно быстрее.
      0
      с замечаниями полностью согласен…
        +3
        Добавьте седьмой код и скажите, на сколько эти советы ускорили — безумно интересно.
          +1
          Можно я попридираюсь? =)
          В шестом листинге закралась ошибка, и в седьмом не исправилась. Появляется еcли n<2.
          lst = n>=2 and [2] or []
            0
            Если уж придираться :)
            lst = [2] if n >= 2 else []
        0
        Только не понятно, почему j * j — 1. Просто j * j. Иначе, если i = 8, мы проверим j = 3, что излишне.
          0
          На самом деле и без int((sqrt(i)+1) можно обойтись. Ведь при последовательном переборе при достаточно больших i значение либо увеличивается на 1, либо остаётся тем же ;)
            0
            Вынос корня за пределы цикла дает бОльшую прибавку к скорости, нежели использование произведения
            В моей реализации это давало ускорение до 10%
            +32
            Поразительно. Вы начали писать программу, не зная про «Решето Эратосфена»?
              0
              про это я знаю и что этот алгоритм работает быстрее тоже в курсе… спасибо… )
                0
                тут можно схитрить, и сказать, что методы на основе решета Эратосферна не могут быть реализованы как генераторы (без ограничения предела генерации) и жрут память при больших n :P
                  +3
                  Сказать можно. А можно подумать как сделать так, чтоб не жрало много памяти. Это я намекаю, что автору стоило именно в этом ключе думать ;)
                  В «Решете Эратосфена» без потери скорости алгоритма можно разбить большой диапазон на куски фиксированной длины. Например, по 1000 элементов. И найдя все простые от 1 до 1000, далее просеивать элементы от 1000 до 2000 и т.д.

                  Хотя об этом статья уже есть на хабре.
                0
                Тоже по вашей теме.
                stevegilham.blogspot.com/2009/01/pipe-operator-in-scala.html

                Не знаю, правда, как по скорости работы. Но выглядит хорошо.
                  0
                  Взято из поста на хабре, ссылку на который потерял.
                  +2
                  хабра скоро станет одной большой лекцией по математике, как например в данном случае по теории чисел )
                    +6
                    Вы так об этом говорите, будто это плохо.
                    Хабр — для It, и алгоритмы — вещь нужная.
                      +9
                      Не вижу, где это я говорю, что это плохо, особенно если учесть, что я по образованию математик
                    +2
                    Готов поспорить, что если убрать проверку на делимость на 5, станет только быстрее. А вот если правильным образом перебирать остатки при делении на 2 * 3 или 2 * 3 * 5, то должно ещё прирасти.
                    И ещё: обычно заводят внутри временную переменную, инициализируюмую как i и при нахождении делителя делят её сколько возможно на этот делитель. Тем самым, корень достигается быстрее.
                      0
                      ∀p > 3 ∈ 6x ± 1
                      0
                      > проверять только те числа, которые заканчиваются на 1, 3, 7 или 9
                      тогда лучше только эти числа и генерить для проверки, для этого можно развернуть цикл и убрать лишнее деление с остатком:

                      def test(i):
                          for j in lst:
                              if j*j-1 > i:
                                  lst.append(i)
                                  break
                              if (i % j == 0):
                                  break
                          else:
                              lst.append(i)
                      
                      m = 1
                      while True:
                          m += 2
                          if m > n: break
                          test(m)
                          m += 4
                          if m > n: break
                          test(m)
                          m += 2
                          if m > n: break
                          test(m)
                          m += 2
                          if m > n: break
                          test(m)
                      
                        0
                        а черт :) там еще 5 должно быть. но идею вы поняли
                          0
                          т.е. случаи n <= 5 можно проверять отдельно
                          +1
                          Я примерно о том же. Но для эффективности будет невредно функцию заинлайнить.
                            0
                            угу, надо бы. я просто не хотел чтоб пост разросся, только идею показать быстро :) даже про 5 забыл…
                          +1
                          Так как мне уже успели насрать в карму, то продолжу здесь писать. Интересны в данном контексте были бы полиномы, которые выдают простые числа. Сразу спешу расстроить: такого, который выдаёт только простые числа, не существует (Гольдбах, 1752). Эйлер же нашёл в 1772 одну совершенно необычную функцию: f_q(x)=x^2+x+q для всех x < q-1 если q из 2, 3, 5, 11, 17, 41.
                            +4
                            Также: любой алгоритм, основанный на сите Эратосфена, перестаёт быть эффективным для больших простых чисел, т.е. для действительно интересных чисел, нужных в криптографии. Для того, чтобы найти алгоритмы для настоящих простых чисел, приходится сначала прорубаться через понятия псевдопростых чисел и чисел Кармихаэля, через тесты для псевдопростых чисел Соловай-Штрасена, Рабин-Миллера итп.
                            Сегодня самыми эффективными можно считать тесты на простые числа при помощи эллиптических кривых, также известны APRCL и AKS. Подробнее писать в комментариях смысла, наверное, нет.
                              +1
                              Вы путаете области применимости алгоритмов: решето Эратосфена ищет все простые числа до некоторой границы, и именно эта задача обсуждается в посте; умные слова типа «псевдопростые числа» относятся к проблеме тестирования заданного числа на простоту.
                              +4
                              Один из интереснейших фактов о простых числах, наверное, этот: между натуральными n и 2n всегда найдётся простое число.

                              Можно сказать так: существуют простые двоичные числа любой длины!

                                +1
                                Но на самом деле всё намного круче: если вы покопаетесь в простых числах, то вы начнёте замечать, что они, ни смотря ни на какую случайность, довольно регулярны, и окажетесь правы конечно же.

                                Чейбышев доказал, что для любого n можно найти такие C_1 и C_2, что
                                C_1*n/log(n) < pi(n) < C_2*n/log(n),
                                где pi(n) это просто количество простых чисел меньше n.

                                Что показывает взаимосвязь натурального логарифма и простых чисел. Однако это ещё не всё! С помощью Чейбышева (но не только) можно доказать следующее:
                                lim_(x -> inf) pi(x)*log(x)/x = 1
                                Называется это Теорема о распределении простых чисел, на неё можно наглядно взглянуть тут: demonstrations.wolfram.com/ThePrimeNumberTheorem/.
                                  +6
                                  Умора. Чейбышев. Это ж надо.
                                  Гражданин математик, прекратите монолог, состоящий из изложения элементарных фактов. Из википедии желающий найдёт всё перечисленное, а так же как пишется фамилия этого учёного.
                                    +1
                                    Sorry, учил все эти вещи не в русской литературе :(
                                      –1
                                      К тому же, по-моему как раз информация для человека, пишущего алгоритмы с ситом Эратосфена, самая подходящая.
                                        0
                                        Я бы не отказался почитать про интересные факты и теоремы из теории чисел, только желательно в виде статьи. В вики это не так просто, когда не знаешь что ищешь :)
                                          +1
                                          Специально зашёл сюда. Там написано всё вышеперечисленное со ссылками на описание.
                                          Вцелом по т.ч. вы, наверное, правы.
                                            0
                                            Черт, сколько гулял по вики по всяким теоремам, а на страничку про простые числа не додумался зайти ) Спасибо
                                            +1
                                            Теория чисел — большая наука, там фактов и теорем много… Про распределение простых чисел с картинками на русском хорошо написано здесь: ega-math.narod.ru/Liv/Zagier.htm.
                                              +2
                                              Лучше рыться не в Вики, а в Mathworld. Например, Number Theory. Материалов и подробностей больше, но тексты сжатые и могут показаться слишком трудными для чтения.

                                              Есть еще The Prime Pages, хотя вы их, наверное, уже видели по ссылкам из Вики. Там изложение достаточно доступное и темы менее специальные, чем на Mathworld.

                                              Ну или скажите, что вам интересно по теории чисел — это моя специальность, могу подготовить материал. Если карма будет. Могу про эффективные алгоритмы некоторых теоретико-числовых вычислений, например. Не очень заумных, честно-честно. Могу про распределение простых чисел в контексте алгоритмов и оценки сложности алгоритмов.
                                                0
                                                Было бы очень классно алгоритм Берлекампа и Ленстра-Ленстра-Ловаса «на пальцах». Или ссылку, где они понятно объясняются. Спасибо!
                                            +1
                                            «Для любого n можно найти такие C_1 и C_2, что ...» — это не нуждается в доказательстве. Чебышёв доказал, что можно найти такие C_1 и C_2, что для любого (достаточно большого) n выполнены неравенства C_1*n/log(n) < pi(n) < C_2*n/log(n) — почувствуйте разницу. Чебышёву не удалось доказать, что предел существует; это сделали уже после него с использованием тяжёлого матана ТФКП.
                                              0
                                              Вы совершенно правы! сам удивлён, что сформулировал теорему настолько неправильно.
                                            0
                                            Сейчас известны намного более сильные факты: например, для достаточно большого n между n и n+n^0.525 всегда найдется простое число.
                                            +3
                                            Формулируйте чётче, в математике, равно как и в программировании, небольшое изменение слов может диаметрально изменить смысл и верность утверждения. Многочлен, выдающий в точности все простые числа, существует, с поправкой на то, что рассматривать нужно только положительные его значения; только такой многочлен зависит от 26 переменных.
                                              0
                                              Да, и этот многочлен не менее красив, чем теоремы Матиясевича. Тоже Сизого читали?)
                                              Это я к тому, что интересующимся теорий чисел не помешает эта книга: Сизый С. В. Лекции по теории чисел.
                                              0
                                              «и тут Остапа понесло...»(с)
                                              +1
                                              А «решето Эратосфена» для чего? Зачем изобретать велосипед?

                                              Не, я конечно понимаю, что всем всегда хочется написать своё. Но зачем, зачем?!
                                                +3
                                                Если бы Эратосфен когда то не хотел придумать что то своё, то сейчас мы бы не знали что существует «решето Эратосфена»…
                                                  +4
                                                  Если мы постоянно будем изобретать велосипед, то 360 км/ч не разгоним никогда.

                                                  Я тоже недоумеваю, зачем этот алгоритм, если есть решето Эратосфена и современное решето Аткина. А тут и объяснено всё, и реализовано.

                                                  А сейчас, чтобы придумать что-то своё, нужно прежде много изучить, а потом уже сделать что-то лучше. Поверьте, гениально простые идеи на низких уровнях потихоньку исчерпались, и сейчас они гениально просты только для людей с достаточным уровнем знаний в той или иной области. Теорему Ферма тоже пыталось доказать много дилетантов (фермисты), толк от этого нулевой.

                                                  зы. Писать алгоритмы на Python — тенденция хорошая. Люблю этот язык.
                                                0
                                                Вот тут в 10 пунтке приводится пример нахождения простых чисел при помощи алгоритма Решето Эратосфен на разных языках программирования.

                                                Вот пример для php:
                                                for($i=range(2,100);($a[]=array_shift($i))&&count($i)>0;)
                                                foreach($i as $k=>$v) if($v%end($a)==0) unset($i[$k]);
                                                
                                                  0
                                                  В тексте для поступления в Computer Science Center была задача на программирование в которой нужно было искать простые числа. Было интересно узнать, как много есть алгоритмов для нахождения их. Но там верхний предел был 109.
                                                    +1
                                                    Sorry за не в тему, тоже хочу поучавствовать cо своим APL:
                                                    image
                                                    Данное решение использует матрицу, содержащую все числа, кроме простых (видимо, то самое решето Эратосфена)
                                                      +1
                                                      Мне кажется, что вместо того, чтобы перебирать только те числа, которые оканчиваются на 1,3,7,9 (а это 2/5 всех чисел) будет еще немного эффективнее перебирать все числа вида 6k(±)1 — все простые числа кроме 2 и 3 имеют эту форму, их просто генерировать и их всего 1/3 от всех чисел.

                                                        0
                                                        А если продолжить эту логику дальше, приходим к решету Эратосфена.
                                                          0
                                                          Нет, это дополняет решето Эратосфена. Смотри мой комментарий ниже
                                                          +1
                                                          В догонку решето Эратосфена на Haskell:

                                                          eratosthenes (x:xs) = x:eratosthenes (filter ((/= 0).(`mod` x)) xs)
                                                          primes = eratosthenes 2 : 3 : [ 6*n+z | n <- [1..], z <- [-1,1] ]
                                                            0
                                                            решето Эратосфена

                                                            Я его на ASM на спектруме как то делал, весёлая штука, но вот с практической точки зрения увы не самый «смак», слишком огромен простор для оптимизации
                                                            0
                                                            Хотя решето Эратосфена позволяет эти числа отсеить, но тут выгода в том, что их можно генерировать (что значительно сократит перебор) и уже затем сгенерированный ряд передавать как основу для решета Эратосфена
                                                            0
                                                            Слово «решето» уже 11 раз встретилось в комментариях, поэтому решил набросать этот алгоритм на Python и сравнить… Листинг 8 добавил в пост… Начиная, примерно с n = 140 «решето» работает быстрее…
                                                              0
                                                              Довольно быстрое решето на питоне, жаль, только, что не генератор.
                                                              def sieve(n):
                                                                      s = xrange(3, n + 1, 2)
                                                                      r = set(s)
                                                                      [r.difference_update(xrange(n << 1, s[-1] + 1, n)) for n in s if n in r]
                                                                      return r.union([2])


                                                              Вот статья на wikibooks про генерацию праймов на разных языках.

                                                              ПС. Кстати, там рядом есть тема, где человек страдает прям как вы.
                                                              0
                                                              Это топик про нахождение простых чисел, так? Тогда добавьте оценку асимптотиечкой сложности к вашим алгоритмам, а питоновые хаки пишите как отдельный топик — они не интересны.
                                                                0
                                                                Спасибо за конструкцию for-else.
                                                                  0
                                                                  Безотносительно алгоритма…
                                                                  Я в свое время (еще в школе) обнаружил такую закономерность — начиная с простого числа 11, путем прибавления 30 получаем простое число. Уже не помню, но в этой закономерности были очень редкие исключения, например 19 + 30 = 49 = 7 * 7.
                                                                    +2
                                                                    Ваша закономерность объясняется тем, что числа вида 30k+p, где p — выбранное простое, большее 5, гарантированно не имеют делителей меньше 7, которые и высевают основную массу (~75%) составных чисел. В таких арифметических прогрессиях (30k+p) всегда бесконечно много простых чисел (теорема Дирихле), хотя при больших k большинство элементов окажутся все же составными.
                                                                    0
                                                                    а почему без yield?
                                                                      –1
                                                                      я не пойму, сегодня день абитуриентов на хабре?
                                                                      первый курс, первый семестр, первый практикум: нахождение НОД, НОК и простых чисел в диапазоне. и это не говоря уже о том, что давно найден эффективный алогритм поиска простых чисел

                                                                      Only users with full accounts can post comments. Log in, please.