Unix Review Column 26 -- Вычисление простых чисел

Randal Schwartz

Перевод Anton Petrusevich <casus@mail.ru> и Alex Ott <ott@phtd.tpu.edu.ru>

Нет никаких сомнений в том, что Perl великолепен в склеивании вместе внешних программ и данных, но Perl также является хорошим средством для проблем чистого программирования, которые вы можете найти в начале курсов по программированию. Например, давайте взглянем на общую задачу вашего первого года обучения, нахождение простых чисел.

Как я вспоминаю из своих ранних занятий математикой, простым числом является любое целое число больше единицы, которое делится только на единицу и само себя, но не на другие целые числа. Предупреждение: я не математик-идиот, так что это определение не является точным, я надеюсь, что кто-то поправит меня.

Так что, давайте выдадим простые числа пользуясь данным определением:

    GUESS: for (my $guess = 2; $guess <= 10000; $guess++) {
      for (my $divisor = 2; $divisor < $guess; $divisor++) {
        next GUESS unless $guess % $divisor;
      }
      print "$guess\n";
    }

Внешний цикл, названный GUESS, нужен для того, чтобы мы могли делать переходы вперед со связью этого цикла, даже изнутри внутреннего цикла. Переменная $guess содержит число, которое мы проверяем ы каждом из циклов. Я определяю для этого цикла лексическую локальную переменную, задавая ей значений равное 2 и увеличивая ее на 1 до тех пор, пока она не достигнет значения 10,000.

Внутренний цикл старается найти все возможные делители начиная с 2 и до того числа, которое мы проверяем. Выражение $guess % $divisor является истинным в том случае, когда имеется остаток от деления, означая, что $guess не делится на $divisor.

На моей машине эта программа выполняется примерно 13 секунд только для того, чтобы получить все простые числа меньшие 10000. Очевидно, что получение 100,000 простых чисел с помощью этой программы будет очень дорогостоящим.

Но давайте применим немного математики к процессу выбора чисел, для того чтобы ускорить его. С одной стороны, глупо пытаться использовать любой делитель, который больше чем квадратный корень искомого числа. (Hаибольший делитель должен иметь другой множитель, который меньше чем квадратный корень, для того, чтобы их произведение было меньше чем искомое число). Так, что давайте немного изменим внутренний цикл:

    GUESS: for (my $guess = 2; $guess <= 10000; $guess++) {
      for (my $divisor = 2; $divisor * $divisor <= $guess; $divisor++) {
        next GUESS unless $guess % $divisor;
      }
      print "$guess\n";
    }

В этом примере заметьте, что было изменено заключительное условие внутреннего цикла. так что мы отбрасываем проверку $divisor, когда делитель имеет значение большее квадратного корня потенциального простого числа, находящегося в переменной $guess.

Hа моей машине этот прием уменьшает время выполнения примерно на пол секунды. Если я увеличу верхний предел до 100000 то я получу выигрыш в 8 секунд. Давайте попробуем выполнить другую оптимизацию кода. После простого числа равного 2 нет простых четных чисел, так что не зачем пытаться использовать четные делители. (Как сказано в моей книге по математике ``доказательство этого оставим читателю''). Так что давайте немного изменим цикл:

    print "2\n";
    GUESS: for (my $guess = 3; $guess <= 100_000; $guess += 2) {
      for (my $divisor = 3; $divisor * $divisor <= $guess; $divisor += 2) {
        next GUESS unless $guess % $divisor;
      }
      print "$guess\n";
    }

Здесь я явно выдаю число 2, поскольку это сделает алгоритм более простым, чем при включении этого числа внутрь цикла. Затем я пробую каждое нечетное число, начиная с 3. Аналогичным образом внутренний цикл пробует использовать в качестве возможного делителя нечетные числа начиная с 3, но все равно останавливается когда квадрат делителя превышает возможное простое число.

И как я предполагал, этот прием уменьшил время выполнения почти вдвое, уменьшившись примерно до 5 секунд. Он все еще не хорошо масштабируется при получении простых чисел меньших миллиона или 10 миллионов. Для того, чтобы сделать это, нам нужно использовать что-то вместо грубой силы.

Следующим крупным исправлением будет введение классического ``Решета Эратосфена'' для нахождения простых чисел. Для того, чтобы сделать это, мы сначала создадим массив, заполненный нулями, размером равным максимальному числу, которое мы хотим найти. Затем, начиная с числа 2, мы присвоим ненулевое значение всем индексам массив кратным 2, показывая, что соответствующее число не является простым. Переходя далее, мы можем выполнить то же самое для числа 3, устанавливая все индексы массива кратные 3 в ненулевое значение. Это выглядит вот так:

    my $UPPER = 100_000;
    my @sieve = (0) x ($UPPER + 1);
    GUESS: for (my $guess = 2; $guess <= $UPPER; $guess++) {
      next GUESS if $sieve[$guess];
      print "$guess\n";
      for (my $mults = $guess * 2; $mults <= $UPPER; $mults += $guess) {
        $sieve[$mults] = 1;
      }
    }

Верхняя граница возможных простых чисел нужны в нескольких местах, так что я поместил это значение в переменную $UPPER. имя этой переменной набрано прописными буквами, что традиционно означает, что мы используем постоянную переменную, установленную как константу для всей программы. Далее, массив @sieve инициализируется нулями начиная с $sieve[0] и до $sieve[$UPPER]. Поскольку размер этого списка на единицу больше значения $UPPER, то мне необходимо добавить 1. Оператор x воспроизводит заданный слева список (просто 0) столько раз, чтобы сделать @sieve достаточно длинным.

Цикл GUESS снова проверяет возможные значения начиная с 2 и до верхнего предела. Если $sieve[$guess] уже установлен, то он уже был установлен предыдущим кандидатом и мы используем умножение целого числа на простое число, найденное раннее, так что мы просто перескакиваем к следующему варианту. В противном случае наш вариант печатается и производится заполнение массива.

Для всех чисел кратных найденному простому числу выполняется проход по массиву, устанавливая индексы с кратными значениями в значение 1. Заметьте, что сдесь нет умножений или делений (если вы не обращаете внимание на $guess * 2, который я мог бы записать как $guess + $guess. Это делает это алгоритм чрезвычайно удобным для процессора.

Как я и предполагал, мы уменьшили время работы до 2 х секунд. При увеличении верхней границы расчетов в этой программе до миллиона, время работы увеличилось до 15 секунд. У меня начался процесс свопинга, поскольку мы создаем в памяти миллион скаляров. В действительности нам нужен только один бит. К счастью, мы можем прямо представить битовый вектор, используя оператор vec:

    my $UPPER = 1_000_000;
    my $sieve = "";
    GUESS: for (my $guess = 2; $guess <= $UPPER; $guess++) {
      next GUESS if vec($sieve,$guess,1);
      print "$guess\n";
      for (my $mults = $guess * 2; $mults <= $UPPER; $mults += $guess) {
        vec($sieve,$mults,1) = 1;
      }
    }

Здесь переменная $UPPER используется для тех же целей, как и в предыдущей программе, но мы теперь используем одиночный скаляр $sieve для хранения группы битов, первоначально равной нулю.

Остаток программы остается неизменным, но мы теперь мы заменили действия с массивом на операции с vec. vec($sieve,$guess,1) имеет истинное значение, когда соответствующий бит в позиции $guess равен 1. И мы можем использовать ту же самую конструкцию в левой части присваивания, для того, чтобы устанавливать значения бита. Достаточно хитро и компактно.

Достаточно странно, но программа все равно выполняется 15 секунд, но свопирование сильно уменьшилось. Оно все равно кажется высоким, но затем я осознал, что я устанавливаю много уже установленных битов. Hам не нужно затрагивать все кратные $guess числа, поскольку все кратные числа меньшие квадрата $guess уже были установлены ранее!

Исправив это, мы получим

    my $UPPER = 1_000_000;
    my $sieve = "";
    GUESS: for (my $guess = 2; $guess <= $UPPER; $guess++) {
      next GUESS if vec($sieve,$guess,1);
      print "$guess\n";
      for (my $mults = $guess * $guess; $mults <= $UPPER; $mults += $guess) {
        vec($sieve,$mults,1) = 1;
      }
    }

что уменьшает время выполнения программы примерно до 13 секунд. Hе такое уж и большое достижение, но каждый бит помогает нам при увеличении объемов. Все таки мы нашли 78,498 простых чисел!

Если вы хотите узнать о простых числах немного больше (возможно больше, чем вы хотите), то я случайно натолкнулся на страницы по адресу http://www.utm.edu/research/primes/, которая является очень хорошим ресурсом.

Так, что вы можете сейчас увидеть, что проблемы первых годов компьютерной науки могут быть легко решены с помощью Perl, показывая, что Perl является хорошим средством и для теоретических и для практических задач. наслаждайтесь.