[Из песочницы] Решето Сундарама
Решето Сундарама в сети представлено большим количеством источников краткой справочной информации. Тем не менее, я решил изложить то, что хотел бы прочитать сам в начале изучения теоретико-числовых алгоритмов.
Решето Сундарама входит в тройку известнейших методов генерации простых чисел. Сейчас к нему принято относиться как к некоторой экзотике по причине плохой вычислительной сложности: 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.
