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

Помечено: 

В этой теме 11 ответов, 2 участника, последнее обновление 12 час., 49 мин. назад.

  • Автор
    Сообщения
  • #5835
    @yurmikh

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

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

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

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

    В основе алгоритма лежит следующий факт. Всякое нечетное составное число n представимо в виде произведения двух натуральных нечетных чисел больше единицы:
    $$ n = (2*i + 1)*(2*j + 1), (1)$$
    где i и j – натуральные числа (ноль не является натуральным числом). Простое число n в таком виде представить невозможно, так как в противном случае это означало бы, что n имеет дополнительные делители кроме самого себя и единицы.

    Запишем нечетное n в виде 2*m + 1, подставим в (1), и получим выражение для m:
    $$ m = 2*i*j + i + j. (2)$$
    Это преобразование непосредственно приводит к идее решета Сундарама.

    Для того чтобы отсеять составные числа из заданного интервала, следует использовать массив, в котором элемент с индексом 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() из конструктора класса. После его окончания поле public double dt будет содержать время выполнения Sieve в секундах.

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

  • #5836
    @admin

    Однако асимптотика – асимптотикой, а на практике в 32-битном диапазоне просеивания Аткин, например, перегоняет Сундарама только при тщательной оптимизации.

    Посмотрел асимптотику решета Аткина — O(N/(log log N)). Это все-таки значительно лучше чем O(N*log(N)). Википедия говорит, что у Аткина лучше оценка по памяти — примерно N^(1/2), когда у Сундарама O(N). Заметить разницу, я думаю, можно будет при N равном 100000.

    • #5876
      @yurmikh

      Все правильно, асимптотика у решета Аткина значительно лучше. Вот только грамотно запрограммировать этот метод не очень просто. Видимо те, кто знает как это сделать, не опускаются до массовой аудитории. Ни разу не видел в сети реализацию реального решета Аткина.

      Тот код, который приводится в Википедии, и варианты, распространенные в сети, просто иллюстрируют возможность использования квадратичных форм, а алгоритм имеет вычислительную сложность O(n) да еще и большой постоянный коэффициент. По этой причине этот алгоритм не выигрывает у решета Сундарама в указанном мной диапазоне.

      Коэффициент можно уменьшить, улучшить время почти в два раза и обогнать Сундарама, но для этого нужно проявить изобретательность. Собирался написать об этом, но от намерения отказался.

      • #5932
        @admin

        Если не секрет, почему вас интересовала эта (Решето Сундарама) тема? — вы где-то применяете такие алгоритмы?

        • #5961
          @yurmikh

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

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

          Удручает беспомощность статей типа этой. На решето Аткина я потратил достаточно много времени, но, к сожалению, не могу написать о нем также ясно, как о решете Сундарама.

          • #5962
            @admin

            По таким темам, скорее всего, стоит искать информацию в научных статьях, а не на хабре.

          • #5963
            @yurmikh

            Реализация — это вне науки. Гугл и Яндекс ничего интересного не находят.

          • #5964
            @admin

            Я тоже решил поискать информацию по решету Аткина :). Посмотрел научные статьи и, в целом, хабр оказался не так уж плох (на их фоне)… Например:

            1) Дмитриев Е. А. в статье «Решето Аткина для поиска простых чисел» из журнала «Приоритетные направления развития образования и науки» копипастит откуда-то алгоритм (как и в остальных своих статьях, кстати). Инетересен график к этой статье:

            Вероятно, что никакой реализации алгоритмов он вообще не делал, а взял формулы асимптотической сложности, вбил их в эксель и… получил графики. \

            2) Березина В. А. «Анализ скорости поиска простых чисел с использованием различных алгоритмов и методов распараллеливания средствами OpenMP» из журнала «Студенческая наука для развития информационного общества» — позорная статья, без реализации и (судя по графикам и выводам) с неправильным распараллеливанием. Сам подход с хаотичной («методом тыка») вставкой директив типа static, dynamic — антинаучный. Не указано на каком железе они это делали. Убило несоответствие данных в таблицах и на графиках. А еще — опция распараллеливания «без параметров», результаты которой магическим образом на 10% отличаются от static (которая и ставится «по умолчанию»). Ну там руководитель (М. В. Трофимова) — канд. пед. наук — так что все нормально.

            3) Сексембаева М.А. в статье «Исследование оптимизации вычислительных систем» (Журнал «Интернаука» No 3 (3). Часть 1, 2016 г. описывает параллельную версию алгоритма Аткина. Сравнивает ее с решетом Сундарама. Ссылки на реализацию там нет, но судя по выводам она смогла реализовать решето Аткина значительно эффективнее. Не исключено (ведь код не видим), что она просто отвратительно реализовала решето Сундарама — поэтому получилось так…

            Таких статей я просмотрел десяток, однако, нашел и что-то более-менее интересное:

            1) Соловьёв В. Г. в статье «Алгоритм решета простых чисел» предлагает новый алгоритм (как я понял — модификация решета Аткина). Мне понравилось описание. Сложилось ощущение, что автор реализовал все это.

            2) В. А. Минаев и др. в статье «Быстродействие индексных алгоритмов при поиске простых чисел по сравнению с решетом Аткина» подробно описывают возможные модификации алгоритма. Опять же, (1) есть ощущение, что они все это реализовали; (2) у этих ребят есть множество других статей по этим темам.

            Я на вашем месте (если вы интересуетесь темой) — написал бы этим ребятам на почту, она есть в статьях.

            Также, я нашел реализации:

            1) В статье «Sieve Of Atkin» (https://fylux.github.io/2017/03/16/Sieve-Of-Atkin/) есть ссылки на реализации на C++ и Maple.

            2) Еще одна реализация: https://github.com/rmoritz/prime-sieve-cxx. К ней тоже есть статья (не научная), но мне она понравилась меньше чем реализация.

            3) Тут: https://github.com/andreasfrom/primekincpp автор честно пишет, что это «его первый опыт» реализации таких алгоритмов, но он пробовал его распараллелить (ИМХО, как-то странно он это делал).

            4) Эта реализация мне понравилась: https://github.com/jinscoe123/atkin — вы более компетентны в вопросе, быть может прокомментируете.

            5) Даже не пытался внивнуть в эти исходники: https://github.com/ESGorelov/Atkin-MoreinPrimeTest/tree/master/Atkin-MoreinPrimeTest , но это чей-то дипломный проект. Наверное, там должна быть «серьезная» работа — не лабораторная, все-таки.

            Другие реализации мне совсем не понравились — они примерно как в википедии.

          • #5975
            @yurmikh

            Первый раз получил сообщение о вашем комментарии на почту. Я обязательно отвечу по существу немного позже. Если вас тема заинтересовала, нужно прочитать статью Аткина здесь. Могу прислать грубый «подстрочник» — перевод. Источником проблем с реализацией является конспективное изложение. Видимо мало кто может понять четвертую часть этой статьи.

          • #5985
            @yurmikh

            Со статьями я знаком, а вот ссылки на реализации для меня оказались новыми. Спасибо. Я просмотрел их «по диагонали». В чужом коде буду ковыряться долго, но подозреваю, что решение своей проблемы не найду. Она заключается в том, что Аткин для поиска решений своих квадратичных форм использует алгоритм Брезенхема. Правильнее, видимо, сказать метод, так как алгоритмов Брезенхема и их модификаций существует много. Метод состоит в проходе по ближайшим к кривой точкам целочисленного растра с целью ее отрисовки. Так вот мне не удается перебросить мостик от Брезенхема к тому, что делает Аткин.

          • #5986
            @yurmikh

            (5) — дипломный проект к просеиванию простых отношения не имеет. Это какой-то самопальный тест на простоту чисел. Я, как и вы, не разбирался.
            Остальные проекты на гитхаб мне знакомы. Я удалил из своей коллекции ссылки на них и поэтому не сразу узнал.
            (1) – это ученическая работа. Автор указывает, что реализация является иллюстрацией. Здесь структурно выделены три алгоритма (из шести) статьи Аткина. Реализация чудовищной неэффективности. Я не буду утверждать, но подозреваю, что некоторые простые будут пропускаться.
            Все остальные основаны на полном переборе значений x и y. (Т.е. для меня интереса не представляют.) Детали я не анализировал.

          • #5987
            @admin

            Спасибо за интересную информацию.

Для ответа в этой теме необходимо авторизоваться.