Без названия и без подзаголовка

Простые числа. Часть первая. Неверный путь.

Простым называется натуральное число, больше единицы, которое делится только на единицу и на само себя. Несложное определение, не так ли? Давным давно Евклид доказал что таких чисел бесконечное множество, так что давайте напишем функцию, которая будет возвращать все простые числа меньше заданного, используя наиболее простой и очевидный метод. В чем он заключается? Нет ничего проще чем делить каждое число на каждое целое большее единицы и меньшее чем само число. Если в процессе деления мы ни разу не получим нулевой остаток, то это будет означать что число простое. Запишем это в Python:

from time import time
def primeList1(n, verbose = True) :
    t = time()
    ret = []
    for j in range(2,n+1) :
        jIsPrime = True
        for k in range(2,j) :
            if j % k == 0 :
            jIsPrime = False
            break
        if jIsPrime :
            ret += [j]
    if verbose :
        print "Calculated primes to",n,"in",time()-t,"sec."

    return ret

Обратите внимание на слабую попытку оптимизации кода — использование инструкции break, которая передаст управление первой операции после цикла, как только будет найден первый делитель. Если вы попробуете выполнить этот код, то он будет вполне сносно работать. Используйте его для получения всех простых чисел меньше 1000 — результат вы получите, скорее всего, мгновенно. Получение всех простых чисел до 10000 затянется на одну-две секунды. Волноваться пока не о чем. Если вы попробуете запустить программу для всех простых меньше 100000, то пока она будет выполняться вы сможете выпить чашечку кофе. А добавив еще несколько нулей, скрипт будет выполняться всю ночь. Это неудивительно. Например, для проверки того что число 9973 простое, нужно проверить все делители от 2 до 9972, включительно. Получается что время выполнения растет как n2.

Есть несколько очень простых усовершенствований, которые могут резко ускорить процесс, и которые требуют немного, совсем немного, математических знаний…

Во-первых, для определения простоты 9973 мы не должны рассматривать все целые из диапазона 2..9972. Мы можем остановить проверку после первого целого, большего чем квадратный корень из 9973. Почему? Потому что если рассматриваемое число делится на целое, большее чем квадратный корень из этого числа, то оно должно делиться и на частное от этого деления, которое очевидно меньше чем квадратный корень. И так как делить на меньшие числа мы уже пытались (безуспешно), то продолжать смысла нет.

Итоговый код будет очень похож на предыдущий:

from time import time
from math import sqrt

def primeList2(n, verbose = True) :
    t = time()
    ret = []
    for j in range(2,n+1) :
        jIsPrime = True
        kMax = int(sqrt(j))
        for k in range(2,kMax+1) :
            if j % k == 0 :
                jIsPrime = False
                break
        if jIsPrime :
            ret += [j]
    if verbose :
        print "Calculated primes to",n,"in",time()-t,"sec."

    return ret

Попробуйте выполнить его, и почувствуйте волнение от новой скорости алгоритма. Если вы попробуете найти все простые числа до 10000, то новая функция выдаст их вам скорее всего… в сто раз быстрее! Если вы будете пробовать ее для больших чисел, то опять получите увеличение скорости, потому что теперь время нашего алгоритма растет со скоростью n1.4.

Мы все еще делаем много лишних проверок. Подумайте опять о числе 9973. Если мы попробуем разделить его на 2 и получим ненулевой остаток, тогда нет смысла проверять делимость на 4, 6, 8 и на все целые, кратные 2. Итак, мы можем удвоить скорость алгоритма, проверяя только нечетные числа, но давайте еще немножко подумаем прежде чем вернемся к кодированию… После неудачи с 2, мы переключимся на 3 и остаток опять будет ненулевой. Тогда какой смысл пробовать 6, 9, 12 и все остальные числа, кратные 3? Если вы продолжите рассуждать в том же духе, то постепенно придете к заключению о том что проверки на делимость нужно делать только для простых чисел, потому как со временем вы начнете проверять делимость 9973 на составные числа. К этому моменту будут проверены все их простые множители, и если прошлые проверки оказались неудачными, то исследуемое число будет простым.

А теперь самое интересное! Разве проверка того является число простым или нет не влияет на то что мы делаем после проверки? Как мы получать простые числа для проверки очередного числа? В самой первой оптимизации мы уже видели что проверочные деления должны продолжаться только до ближайшего целого, которое больше корня из проверяемого числа. Получается что нам нужно знать только первое простое число — 2. Далее мы можем использовать массив уже проверенных простых чисел для формирования списка элементов для проверки… Код, который будет отвечать всему вышесказанному будет выглядеть примерно так:

from time import time
from math import sqrt

def primeList3(n, verbose = True) :
    t = time()
    ret = []
    for j in range(2,n+1) :
        jIsPrime = True
        kMax = sqrt(j)
        for k in ret :
            if k > kMax :
                break
            if j % k == 0 :
                jIsPrime = False
                break
        if jIsPrime :
            ret += [j]

    if verbose :
        print "Calculated primes to",n,"in",time()-t,"sec."

    return ret

Если вы протестируете этот код, то увидите что скорость возросла. Последняя оптимизация не сделала его в сто раз быстрее, но на моем компьютере новый код выполняется в 3,5 раза быстрее чем предыдущий при вычислении всех простых чисел до 100000 и в 5,5 раз быстрее при вычислении всех простых до 100000. Обратите внимание — эффективным код делает не четырехкратное увеличение в скорости, а тот факт что чем больше число на входе, тем лучше работает новый код. Похоже, что теперь время выполнения растет как n1.25

Последняя версия содержит основополагающие идеи по-настоящему «быстрого» кода. Я имею ввиду пропуск произведений простых чисел. Но реализация все еще опирается на медленную и трудоемкую операцию проверки делимости. Избавление от нее — предмет следующей оптимизации, которую вы увидите в следующем посте.

Оригинал: http://numericalrecipes.wordpress.com/2009/03/05/prime-numbers-1-the-wrong-way/