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

Решето Сундарама входит в тройку известнейших методов генерации простых чисел. Сейчас к нему принято относиться как к некоторой экзотике по причине плохой вычислительной сложности: 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 = (v(2*mmax+1)-1)/2 = (vN-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.