Решето Сундарама

Решето Сундарама в сети представлено большим количеством источников краткой справочной информации. Тем не менее, я решил изложить то, что хотел бы прочитать сам в начале изучения теоретико-числовых алгоритмов.

Решето Сундарама входит в тройку известнейших методов генерации простых чисел. Сейчас к нему принято относиться как к некоторой экзотике по причине плохой вычислительной сложности: O(N(logN)). Однако асимптотика – асимптотикой, а на практике в 32-битном диапазоне просеивания Аткин, например, перегоняет Сундарама только при тщательной оптимизации.

Реализации решета Аткина, имеющие хождение в интернете, не превосходят решето Сундарама ни по временным характеристикам, ни по эффективности использования памяти. Так что метод Сундарама вполне можно использовать как вспомогательный инструмент при экспериментах с более продвинутыми алгоритмами.
Решето Сундарама находит все простые числа в некотором заданном диапазоне натуральных чисел 3 ? n ? N «отсеивая» составные. Без потери общности, значение N можно считать нечетным. Если N четно, то оно гарантированно составное, и его можно исключить из диапазона просеивания, уменьшив значение верхней границы на единицу.

В основе алгоритма лежит следующий факт. Всякое нечетное составное число n представимо в виде произведения двух натуральных нечетных чисел больше единицы:

(1) n = (2*i + 1)*(2*j + 1),

где i и j – натуральные числа (ноль не является натуральным числом). Простое число n в таком виде представить невозможно, так как в противном случае это означало бы, что n имеет дополнительные делители кроме самого себя и единицы.

Запишем нечетное n в виде 2*m + 1, подставим в (1), и получим выражение для m:

(2) m = 2*i*j + i + j.

Это преобразование непосредственно приводит к идее решета Сундарама.

Для того чтобы отсеять составные числа из заданного интервала, следует использовать массив, в котором элемент с индексом m является представителем числа 2*m + 1. Организовав перебор значений переменных i и j, будем вычислять индекс m по формуле (2), и в соответствующих элементах массива устанавливать отметку «составное число». После завершения перебора все составные числа диапазона будут отмечены, и простые числа можно выбрать по признаку отсутствия отметки.

Остается уточнить технические детали.

Исходя из предыдущих рассуждений, для представления верхней (нечетной) границы диапазона просеивания N, индекс m принимает свое максимальное значение mmax = (N – 1)/2. Этим определяется необходимый размер массива.

Перебор значений i и j будем осуществлять в двух циклах: внешний цикл для перебора значений i, и вложенный внутренний цикл для значений j.

Начальным значением счетчика внешнего цикла является i = 1. Учитывая полную симметричность вхождения переменных i и j в выражение (2), для исключения повторных дублирующих вычислений, внутренний цикл следует начинать со значения j = i.

Найдем максимальное значение для счетчика внешнего цикла imax ? i. Если граница диапазона N является нечетным составным числом, то при значении i = imax, внутренний цикл должен выполниться как минимум один раз со значением своего параметра j = imax для вычеркивания N, а выражение (2) при этом примет свое максимальное значение:

mmax = 2*imax*imax + imax + imax,
imax^2 + imax — mmax/2 = 0.

Решив это квадратное уравнение, получим: imax = (?(2*mmax+1)-1)/2 = (?N-1)/2.
Ограничение для счетчика внутреннего цикла jmax ? j найдем из неравенства
mmax ? 2*i*j + i + j, откуда получаем: jmax = (mmax – i)/(2*i + 1).

Значения верхних пределов следует округлять до целого значения, отбрасывая дробную часть.

Ниже показан листинг очень простого класса Sundaram на языке C#, реализующий описанный метод.

using System; using System.Collections;  namespace CSh_Sundaram {     public class Sundaram     {         public double dt;          // время просеивания (сек)         private long t1;           // начало просеивания         private long t2;           // окончание просеивания         private uint limit;        // верхняя граница диапазона просеивания         private int arrlength;     // размер массива         private BitArray Prim;     // массив результатов просеивания         private int counter;          public Sundaram(uint _limit)         {             limit = _limit;             if (limit % 2 == 0) limit -= 1;             arrlength = (int)(limit / 2);             Prim = new BitArray(arrlength);             t1 = DateTime.Now.Ticks;             Sieve();                 // Просеивание             t2 = DateTime.Now.Ticks;             dt = (double)(t2 - t1) / 10000000D;             counter = -1;         }          public uint NextPrime         {             get {                 while (counter < arrlength - 1) {                     counter++;                     if (!Prim[counter])                         return (uint)(2 * counter + 3);                 }                 return 0;             }         }          private void Sieve() {             int imax = ((int)Math.Sqrt(limit) - 1) / 2;             for (int i = 1; i <= imax; i++) {                 int jmax = (arrlength - i) / (2 * i + 1);                 for (int j = i; j <= jmax; j++) {                     Prim[2 * i * j + i + j - 1] = true;                 }             }         }     } } 

В качестве параметра при создании объекта типа Sundaram указывается верхняя граница диапазона просеивания. Для просеивания класс использует массив типа BitArray. Это увеличивает время выполнения, но позволяет перекрыть весь 32-разрядный диапазон типа uint. Просеивание выполняется при обращении к методу Sieve() из конструктора класса. После его окончания поле dt будет содержать время выполнения Sieve в секундах.

При обращениях к свойству NextPrime последовательно, начиная со значения 3, в порядке возрастания, будут возвращаться простые числа. После того, как все простые из диапазона исчерпаны, будет возвращаться значение 0.

Источник: habr.com

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