C #: Как сделать Сито Аткина инкрементным

Я не знаю, возможно ли это или нет, но я просто должен спросить. Мои математические и алгоритмические навыки меня здесь не устраивают: P

Дело в том, что у меня теперь есть этот class, который генерирует простые числа до определенного предела:

public class Atkin : IEnumerable { private readonly List primes; private readonly ulong limit; public Atkin(ulong limit) { this.limit = limit; primes = new List(); } private void FindPrimes() { var isPrime = new bool[limit + 1]; var sqrt = Math.Sqrt(limit); for (ulong x = 1; x <= sqrt; x++) for (ulong y = 1; y <= sqrt; y++) { var n = 4*x*x + y*y; if (n <= limit && (n % 12 == 1 || n % 12 == 5)) isPrime[n] ^= true; n = 3*x*x + y*y; if (n  y && n <= limit && n % 12 == 11) isPrime[n] ^= true; } for (ulong n = 5; n <= sqrt; n++) if (isPrime[n]) { var s = n * n; for (ulong k = s; k <= limit; k += s) isPrime[k] = false; } primes.Add(2); primes.Add(3); for (ulong n = 5; n <= limit; n++) if (isPrime[n]) primes.Add(n); } public IEnumerator GetEnumerator() { if (!primes.Any()) FindPrimes(); foreach (var p in primes) yield return p; } IEnumerator IEnumerable.GetEnumerator() { return GetEnumerator(); } } 

Теперь я хотел бы избавиться от предела, чтобы последовательность никогда не прекращалась (теоретически).

Я думаю, что это может произойти примерно так:

  1. Начнем с некоторого тривиального предела
    • Найти все простые числа до предела
    • Допускайте все новообретенные простые числа
    • Увеличьте предел (удвоение или возведение в квадрат старого предела или что-то в этом роде)
    • Перейти к шагу 2

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

Есть ли способ сделать это? Моя основная проблема заключается в том, что я не совсем понимаю, что такое x и y например, в этом алгоритме. Например, могу ли я просто использовать один и тот же тип алгоритма, но установить x и y в oldLimit (изначально 1) и запустить его до newLimit ? Или как это будет работать? Любые яркие умы с некоторым освещением, чтобы избавиться от этого?

Дело в том, что мне не придется устанавливать этот предел. Так что я могу, например, использовать Linq и просто Take() но многие простые числа, которые мне нужны, не беспокоясь о том, является ли предел достаточно высоким и т. Д.

Вот решение неограниченного простого просеивания в C #, которое может быть реализовано с использованием алгоритмов Sieve of Eratosthenes (SoE) или Sieve of Atkin (SoA); однако я утверждаю, что вряд ли стоит сложная задача оптимизированного решения SoA, чем истинное SoE дает одинаковую производительность без такой сложности. Таким образом, это, пожалуй, только частичный ответ, в то время как он показывает, как реализовать лучший алгоритм SoA и показывает, как реализовать неограниченную последовательность с помощью SoE, он только намекает на то, как их объединить для написания разумно эффективного SoA.

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

Сначала мы должны прокомментировать точку этого упражнения при создании неограниченной последовательности простых чисел, чтобы разрешить использование методов расширения IEnumerable, таких как Take (), TakeWhile (), Where (), Count () и т. Д., Поскольку эти методы отдают некоторые из-за добавленных уровней вызова метода, добавляя по меньшей мере 28 машинных циклов для каждого вызова и возвращаемого значения для enums одного значения и добавления нескольких уровней вызовов методов на каждую функцию; при этом, имея (фактически) бесконечную последовательность, по-прежнему полезно, даже если для получения более высокой скорости использовать более императивные методы фильтрации для результатов перечислений.

Обратите внимание, что в более простых примерах я ограничил диапазон значений беззнаковых 32-разрядных чисел (uint) в прошлом, поскольку базовые реализации SoE или SoA не очень подходят для эффективности, и нужно переключиться к «ковшу» или другой форме инкрементного сита для части просеивания для эффективности; тем не менее, для полностью оптимизированного сеточного сечения страницы, как и в самой быстрой реализации, диапазон увеличивается до 64-битного диапазона, хотя, как написано, вероятно, он не захочет просеять около ста триллионов (от десяти до четырнадцатой мощности) в качестве времени выполнения потребуется до сотен лет для полного диапазона.

Поскольку вопрос выбирает SoA, возможно, неправильные причины вначале ошибочно принимают предварительное сито Trial Division (TD) для истинного SoE, а затем используют наивный SoA над ситом TD, давайте создадим истинное ограниченное SoE, которое может быть реализовано в нескольких строки как метод (который может быть преобразован в class согласно стилю кодирования реализации вопроса) следующим образом:

 static IEnumerable primesSoE(uint top_number) { if (top_number < 2u) yield break; yield return 2u; if (top_number < 3u) yield break; var BFLMT = (top_number - 3u) / 2u; var SQRTLMT = ((uint)(Math.Sqrt((double)top_number)) - 3u) / 2u; var buf = new BitArray((int)BFLMT + 1,true); for (var i = 0u; i <= BFLMT; ++i) if (buf[(int)i]) { var p = 3u + i + i; if (i <= SQRTLMT) { for (var j = (p * p - 3u) / 2u; j <= BFLMT; j += p) buf[(int)j] = false; } yield return p; } } 

Это вычисляет простые числа до 2 миллионов в около 16 миллисекунд на i7-2700K (3,5 ГГц) и 203 280 221 простых числах до 4 294 967 291 в 32-битном диапазоне чисел примерно за 67 секунд на одном компьютере (с учетом запасных 256 мегабайт ОЗУ Память).

Теперь использование версии выше для сравнения с SoA вряд ли справедливо, так как истинный SoA автоматически игнорирует кратные 2, 3 и 5, поэтому введение факторизации колес, чтобы сделать то же самое для SoE, дает следующий код:

 static IEnumerable primesSoE(uint top_number) { if (top_number < 2u) yield break; yield return 2u; if (top_number < 3u) yield break; yield return 3u; if (top_number < 5u) yield break; yield return 5u; if (top_number < 7u) yield break; var BFLMT = (top_number - 7u) / 2u; var SQRTLMT = ((uint)(Math.Sqrt((double)top_number)) - 7u) / 2u; var buf = new BitArray((int)BFLMT + 1,true); byte[] WHLPTRN = { 2, 1, 2, 1, 2, 3, 1, 3 }; for (uint i = 0u, w = 0u; i <= BFLMT; i += WHLPTRN[w], w = (w < 7u) ? ++w : 0u) if (buf[(int)i]) { var p = 7u + i + i; if (i <= SQRTLMT) { var pX2 = p + p; uint[] pa = { p, pX2, pX2 + p }; for (uint j = (p * p - 7u) / 2u, m = w; j <= BFLMT; j += pa[WHLPTRN[m] - 1u], m = (m < 7u) ? ++m : 0u) buf[(int)j] = false; } yield return p; } } 

Вышеприведенный код вычисляет простые числа до 2 миллионов в течение примерно 10 миллисекунд, а числа до 32-разрядного диапазона - примерно 40 секунд на том же компьютере, что и выше.

Далее, давайте определим, действительно ли версия SoA, которую мы, вероятно, будем писать здесь, имеет какую-либо выгоду по сравнению с SoE в соответствии с приведенным выше кодом, насколько скорость выполнения. Следующий код следует за моделью SoE выше и оптимизирует псевдокод SoA из статьи Википедии в отношении настройки диапазонов переменных «x» и «y» с использованием отдельных циклов для отдельных квадратичных решений, как это предлагается в этой статье, решение квадратичных уравнений (и квадратичных свободных исключений) только для нечетных решений, сочетающих квадратичную форму «3 * x ^ 2» для решения как для положительного, так и для отрицательного вторых членов в одном и том же проходе и исключения вычислительно дорогостоящих модульных операций, чтобы производят код, который более чем в три раза быстрее, чем наивная версия псевдокода, размещенного там, и используется здесь в этом вопросе. Код ограниченного SoA следующий:

 static IEnumerable primesSoA(uint top_number) { if (top_number < 2u) yield break; yield return 2u; if (top_number < 3u) yield break; yield return 3u; if (top_number < 5u) yield break; var BFLMT = (top_number - 3u) / 2u; var buf = new BitArray((int)BFLMT + 1, false); var SQRT = (uint)(Math.Sqrt((double)top_number)); var SQRTLMT = (SQRT - 3u) / 2u; for (uint x = 1u, s = 1u, mdx12 = 5u, dmdx12 = 0u; s <= BFLMT; ++x, s += ((x << 1) - 1u) << 1) { for (uint y = 1u, n = s, md12 = mdx12, dmd12 = 8u; n <= BFLMT; y += 2, n += (y - 1u) << 1) { if ((md12 == 1) || (md12 == 5)) buf[(int)n] = buf[(int)n] ^ true; md12 += dmd12; if (md12 >= 12) md12 -= 12; dmd12 += 8u; if (dmd12 >= 12u) dmd12 -= 12u; } mdx12 += dmdx12; if (mdx12 >= 12u) mdx12 -= 12u; dmdx12 += 8u; if (dmdx12 >= 12u) dmdx12 -= 12u; } for (uint x = 1u, s = 0u, mdx12 = 3u, dmdx12 = 8u; s <= BFLMT; ++x, s += x << 1) { int y = 1 - (int)x, n = (int)s, md12 = (int)mdx12, dmd12 = ((-y - 1) << 2) % 12; for (; (y < 0) && (uint)n <= BFLMT; y += 2, n += (-y + 1) << 1) { if (md12 == 11) buf[(int)n] = buf[(int)n] ^ true; md12 += dmd12; if (md12 >= 12) md12 -= 12; dmd12 += 4; if (dmd12 >= 12) dmd12 -= 12; } if (y < 1) { y = 2; n += 2; md12 += 4; dmd12 = 0; } else { n += 1; md12 += 2; dmd12 = 8; } if (md12 >= 12) md12 -= 12; for (; (uint)n <= BFLMT; y += 2, n += (y - 1) << 1) { if (md12 == 7) buf[(int)n] = buf[(int)n] ^ true; md12 += dmd12; if (md12 >= 12) md12 -= 12; dmd12 += 8; if (dmd12 >= 12) dmd12 -= 12; } mdx12 += dmdx12; if (mdx12 >= 12) mdx12 -= 12; dmdx12 += 4; if (dmdx12 >= 12) dmdx12 -= 12; } for (var i = 0u; i<=BFLMT; ++i) if (buf[(int)i]) { var p = 3u+i+i; if (i<=SQRTLMT) { var sqr = p*p; for (var j = (sqr - 3ul) / 2ul; j <= BFLMT; j += sqr) buf[(int)j] = false; } yield return p; } } 

Это все еще более чем в два раза медленнее, чем алгоритм SoE для факсимизации колес, выведенный из-за не полностью оптимизированного количества операций. Дальнейшая оптимизация может быть выполнена в коде SoA, как при использовании операций по модулю 60, так и в отношении исходного алгоритма (не псевдокода) и использования бит-бит для обработки только кратных, исключая 3 и 5, чтобы получить код, достаточно близкий по производительности к SoE или даже немного превзойти его, но мы все ближе и ближе усложняем реализацию Berstein для достижения этой производительности. Моя точка зрения: «Почему SoA, когда очень много работает, чтобы добиться такой же производительности, как и SoE?».

Теперь для последовательности неограниченных простых чисел самый простой способ сделать это - определить const top_number из Uint32.MaxValue и исключить аргумент в методах primesSoE или primesSoA. Это несколько неэффективно в том, что отбраковка все еще выполняется по всему диапазону чисел, независимо от того, обрабатывается фактическое основное значение, что делает определение для небольших диапазонов простых чисел намного длиннее, чем должно было бы. Есть и другие причины перейти к сегментированной версии симулята простых чисел, кроме этого и использования экстремальной памяти. Во-первых, алгоритмы, которые используют данные, которые в основном находятся в кэшах данных L1 или L2 процессора, ускоряются из-за более эффективного доступа к памяти и во-вторых, поскольку сегментация позволяет легко разбивать задание на страницы, которые можно отбирать в фоновом режиме с использованием многоядерных процессоров для ускорения, которые могут быть пропорциональны количеству используемых сердечников.

Для эффективности мы должны выбрать размер массива размера кэша CPU L1 или L2, который составляет не менее 16 килобайт для большинства современных процессоров, чтобы избежать обхода кэша, поскольку мы отбираем композиты простых чисел из массива, а это значит, что BitArray может иметь размером восемь раз больше (8 бит на байт) или 128 килобитов. Поскольку нам нужно только просеивать нечетные простые числа, это представляет собой диапазон чисел более четверти миллиона, а это означает, что сегментированная версия будет использовать только восемь сегментов для решеток до 2 миллионов, как того требует проблема Эйлера 10. Можно было бы сохранить результаты из последний сегмент и продолжить с этой точки, но это не позволит адаптировать этот код к многоядерному случаю, когда один из них отбрасывает фон на несколько streamов, чтобы в полной мере использовать многоядерные процессоры. Код C # для неограниченного SoE (одиночного streamа) выглядит следующим образом:

 static IEnumerable primesSoE() { yield return 2u; yield return 3u; yield return 5u; const uint L1CACHEPOW = 14u + 3u, L1CACHESZ = (1u << (int)L1CACHEPOW); //for 16K in bits... const uint BUFSZ = L1CACHESZ / 15u * 15u; //an even number of wheel rotations var buf = new BitArray((int)BUFSZ); const uint MAXNDX = (uint.MaxValue - 7u) / 2u; //need maximum for number range var SQRTNDX = ((uint)Math.Sqrt(uint.MaxValue) - 7u) / 2u; byte[] WHLPTRN = { 2, 1, 2, 1, 2, 3, 1, 3 }; //the 2,3,5 factorial wheel, (sum) 15 elements long byte[] WHLPOS = { 0, 2, 3, 5, 6, 8, 11, 12 }; //get wheel position from index byte[] WHLNDX = { 0, 0, 1, 2, 2, 3, 4, 4, 5, 5, 5, 6, 7, 7, 7, //get index from position 0, 0, 1, 2, 2, 3, 4, 4, 5, 5, 5, 6, 7 }; //allow for overflow byte[] WHLRNDUP = { 0, 2, 2, 3, 5, 5, 6, 8, 8, 11, 11, 11, 12, 15, //allow for overflow... 15, 15, 17, 17, 18, 20, 20, 21, 23, 23, 26, 26, 26, 27 }; uint BPLMT = (ushort.MaxValue - 7u) / 2u; var bpbuf = new BitArray((int)BPLMT + 1, true); for (var i = 0; i <= 124; ++i) if (bpbuf[i]) { var p = 7 + i + i; //initialize baseprimes array for (var j = (p * p - 7) / 2; j <= BPLMT; j += p) bpbuf[j] = false; } var pa = new uint[3]; for (uint i = 0u, w = 0, si = 0; i <= MAXNDX; i += WHLPTRN[w], si += WHLPTRN[w], si = (si >= BUFSZ) ? 0u : si, w = (w < 7u) ? ++w : 0u) { if (si == 0) { buf.SetAll(true); for (uint j = 0u, bw = 0u; j <= BPLMT; j += WHLPTRN[bw], bw = (bw < 7u) ? ++bw : 0u) if (bpbuf[(int)j]) { var p = 7u+j+j; var pX2=p+p; var k = p * (j + 3u) + j; if (k >= i + BUFSZ) break; pa[0] = p; pa[1] = pX2; pa[2] = pX2 + p; var sw = bw; if (k < i) { k = (i - k) % (15u * p); if (k != 0) { var os = WHLPOS[bw]; sw = os + ((k + p - 1u) / p); sw = WHLRNDUP[sw]; k = (sw - os) * p - k; sw = WHLNDX[sw]; } } else k -= i; for (; k 

Вышеприведенный код занимает около 16 миллисекунд, чтобы просеять простые числа до двух миллионов и около 30 секунд, чтобы пролить полный 32-разрядный диапазон чисел. Этот код быстрее, чем аналогичная версия, используя факториальное колесо без сегментации для больших диапазонов чисел, потому что, несмотря на то, что код более сложный по отношению к вычислениям, очень много времени сохраняется, не разбирая кэширование CPU. Большая часть сложности заключается в вычислении новых начальных индексов на базовый базис на каждый сегмент, который можно было бы сократить, сохранив состояние каждого в расчете на каждый сегмент, но приведенный выше код можно легко преобразовать, чтобы запустить процессы отбора отдельные streamи фона для дальнейшего четырехкратного ускорения для машины с четырьмя ядрами, еще больше с восемью ядрами. Не использование classа BitArray и обращение к отдельным местам бит с помощью бит-масок обеспечивали бы дополнительную скорость примерно в два раза, а дальнейшее увеличение факториального колеса обеспечило бы еще одно сокращение времени примерно до двух третей. Лучшая упаковка бит-массива в исключенных индексах для значений, исключенных с помощью факторизации колес, немного повысит эффективность для больших диапазонов, но также добавит немного сложности в манипуляции бит. Тем не менее, все эти оптимизации SoE не будут соответствовать экстремальной сложности реализации Berstein SoA, но будут работать с одинаковой скоростью.

Чтобы преобразовать вышеприведенный код из SoE в SoA, нам нужно только изменить код отсечения SoA из ограниченного случая, но с модификацией, что начальные адреса пересчитываются для каждого сегмента страницы, как будто начальный адрес вычисляется для случая SoE , но даже несколько более сложными в работе, так как квадратики продвигаются квадратами чисел, а не просто кратными штрихами. Я оставлю необходимые адаптации ученику, хотя я действительно не вижу смысла в том, что SoA с разумной оптимизацией не быстрее текущей версии SoE и значительно сложнее;)

EDIT_ADD:

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

Наиболее эффективная и простая оптимизация - это отключение операций отбраковки на сегментную страницу до фоновых streamов, так что, учитывая достаточное количество ядер процессора, основной предел enums простых чисел - это сам цикл enums, и в этом случае все простые числа в 32-битном номере диапазон может быть перечислит примерно через десять секунд на вышеупомянутом эталонном компьютере (на C #) без других оптимизаций, причем все остальные оптимизации include в себя запись реализаций интерфейса IEnumerable, а не использование встроенных iteratorов, уменьшающих это до примерно шести секунд для всех 203 280 221 простых чисел в 32-битном диапазоне номеров (от 1,65 секунды до одного миллиарда), снова с большей частью времени, затрачиваемым только на перечисление результатов. Следующий код содержит многие из этих изменений, в том числе использование фоновых задач из Dotnet Framework 4 ThreadPool, используемого задачами, с использованием конечного автомата, реализованного в виде таблицы поиска, для реализации дополнительной битовой упаковки, чтобы ускорить перечисление простых чисел и переписать алгоритм как перечислимый class с использованием iteratorов «roll-your-own» для повышения эффективности:

 class fastprimesSoE : IEnumerable, IEnumerable { struct procspc { public Task tsk; public uint[] buf; } struct wst { public byte msk; public byte mlt; public byte xtr; public byte nxt; } static readonly uint NUMPROCS = (uint)Environment.ProcessorCount + 1u; const uint CHNKSZ = 1u; const uint L1CACHEPOW = 14u, L1CACHESZ = (1u << (int)L1CACHEPOW), PGSZ = L1CACHESZ >> 2; //for 16K in bytes... const uint BUFSZ = CHNKSZ * PGSZ; //number of uints even number of caches in chunk const uint BUFSZBTS = 15u * BUFSZ << 2; //even in wheel rotations and uints (and chunks) static readonly byte[] WHLPTRN = { 2, 1, 2, 1, 2, 3, 1, 3 }; //the 2,3,5 factorial wheel, (sum) 15 elements long static readonly byte[] WHLPOS = { 0, 2, 3, 5, 6, 8, 11, 12 }; //get wheel position from index static readonly byte[] WHLNDX = { 0, 1, 1, 2, 3, 3, 4, 5, 5, 6, 6, 6, 7, 0, 0, 0 }; //get index from position static readonly byte[] WHLRNDUP = { 0, 2, 2, 3, 5, 5, 6, 8, 8, 11, 11, 11, 12, 15, 15, 15, //allow for overflow... 17, 17, 18, 20, 20, 21, 23, 23, 26, 26, 26, 27, 30, 30, 30 }; //round multiplier up const uint BPLMT = (ushort.MaxValue - 7u) / 2u; const uint BPSZ = BPLMT / 60u + 1u; static readonly uint[] bpbuf = new uint[BPSZ]; static readonly wst[] WHLST = new wst[64]; static void cullpg(uint i, uint[] b, int strt, int cnt) { Array.Clear(b, strt, cnt); for (uint j = 0u, wp = 0, bw = 0x31321212u, bi = 0u, v = 0xc0881000u, bm = 1u; j <= BPLMT; j += bw & 0xF, wp = wp < 12 ? wp += bw & 0xF : 0, bw = (bw > 3u) ? bw >>= 4 : 0x31321212u) { var p = 7u + j + j; var k = p * (j + 3u) + j; if (k >= (i + (uint)cnt * 60u)) break; if ((v & bm) == 0u) { if (k < i) { k = (i - k) % (15u * p); if (k != 0) { var sw = wp + ((k + p - 1u) / p); sw = WHLRNDUP[sw]; k = (sw - wp) * p - k; } } else k -= i; var pd = p / 15; for (uint l = k / 15u + (uint)strt * 4u, lw = ((uint)WHLNDX[wp] << 3) + WHLNDX[k % 15u]; l < (uint)(strt + cnt) * 4u; ) { var st = WHLST[lw]; b[l >> 2] |= (uint)st.msk << (int)((l & 3) << 3); l += st.mlt * pd + st.xtr; lw = st.nxt; } } if ((bm <<= 1) == 0u) { v = bpbuf[++bi]; bm = 1u; } } } static fastprimesSoE() { for (var x = 0; x < 8; ++x) { var p = 7 + 2 * WHLPOS[x]; var pr = p % 15; for (int y = 0, i = ((p * p - 7) / 2); y < 8; ++y) { var m = WHLPTRN[(x + y) % 8]; i %= 15; var n = WHLNDX[i]; i += m * pr; WHLST[x * 8 + n] = new wst { msk = (byte)(1 << n), mlt = m, xtr = (byte)(i / 15), nxt = (byte)(8 * x + WHLNDX[i % 15]) }; } } cullpg(0u, bpbuf, 0, bpbuf.Length); } //init baseprimes class nmrtr : IEnumerator, IEnumerator, IDisposable { procspc[] ps = new procspc[NUMPROCS]; uint[] buf; Task dlycullpg(uint i, uint[] buf) { return Task.Factory.StartNew(() => { for (var c = 0u; c < CHNKSZ; ++c) cullpg(i + c * PGSZ * 60, buf, (int)(c * PGSZ), (int)PGSZ); }); } public nmrtr() { for (var i = 0u; i < NUMPROCS; ++i) ps[i] = new procspc { buf = new uint[BUFSZ] }; for (var i = 1u; i < NUMPROCS; ++i) { ps[i].tsk = dlycullpg((i - 1u) * BUFSZBTS, ps[i].buf); } buf = ps[0].buf; } public uint Current { get { return this._curr; } } object IEnumerator.Current { get { return this._curr; } } uint _curr; int b = -4; uint i = 0, w = 0; uint v, msk = 0; public bool MoveNext() { if (b < 0) if (b == -1) { _curr = 7; b += (int)BUFSZ; } else { if (b++ == -4) this._curr = 2u; else this._curr = 7u + ((uint)b << 1); return true; } do { i += w & 0xF; if ((w >>= 4) == 0) w = 0x31321212u; if ((this.msk <<= 1) == 0) { if (++b >= BUFSZ) { b = 0; for (var prc = 0; prc < NUMPROCS - 1; ++prc) ps[prc] = ps[prc + 1]; ps[NUMPROCS - 1u].buf = buf; var low = i + (NUMPROCS - 1u) * BUFSZBTS; ps[NUMPROCS - 1u].tsk = dlycullpg(i + (NUMPROCS - 1u) * BUFSZBTS, buf); ps[0].tsk.Wait(); buf = ps[0].buf; } v = buf[b]; this.msk = 1; } } while ((v & msk) != 0u); if (_curr > (_curr = 7u + i + i)) return false; else return true; } public void Reset() { throw new Exception("Primes enumeration reset not implemented!!!"); } public void Dispose() { } } public IEnumerator GetEnumerator() { return new nmrtr(); } IEnumerator IEnumerable.GetEnumerator() { return new nmrtr(); } } 

Обратите внимание, что этот код не быстрее последней версии для небольших диапазонов простых чисел, как до одного или двух миллионов из-за накладных расходов на настройку и инициализацию массивов, но значительно быстрее для более крупных диапазонов до четырех миллиардов плюс. Это примерно в 8 раз быстрее, чем реализация ситча Аткина, но идет от 20 до 50 раз быстрее, чем для более крупных диапазонов до четырех миллиардов плюс. Определены константы в коде, устанавливающие размер базового кеша и сколько из них отбираются на stream (CHNKSZ), который может слегка улучшить производительность. Можно было бы попытаться провести некоторые небольшие оптимизации, которые могли бы сократить время выполнения для больших простых чисел до двух раз, таких как дополнительная упаковка бит, как для колеса 2,3,5,7, а не только 2,3,5 колеса для сокращение примерно на 25% в упаковке (возможно, сокращение времени выполнения до двух третей) и предварительное отбраковка сегментов страницы колесиком факториала до коэффициента 17 для дальнейшего уменьшения примерно на 20%, но это примерно так о том, что можно сделать в чистом коде C # по сравнению с C-кодом, который может воспользоваться преимуществами уникальной оптимизации собственного кода.

Во всяком случае, вряд ли стоит дальнейших оптимизаций, если намереваться использовать интерфейс IEnumberable для вывода, поскольку вопрос требует, так как около двух третей времени используется только для enums найденных простых чисел, и только около трети времени тратится на отбраковку составные числа. Лучшим подходом было бы написать методы classа для реализации желаемых результатов, таких как CountTo, ElementAt, SumTo и т. Д., Чтобы полностью исключить перечисление.

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

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

  1. Диапазон использования был увеличен до 64-но неподписанного числа в диапазоне 18 446 744 073 709 551 615 с удалением проверок переполнения диапазона, поскольку маловероятно, что кто-то захочет запустить программу в течение сотен лет, которые потребуется для обработки всего диапазона чисел до это ограничение. Это связано с очень небольшими затратами времени обработки, поскольку пейджинг может быть выполнен с использованием 32-битных диапазонов страниц, и только конечный первичный вывод должен быть вычислен как 64-разрядное число.
  2. Он увеличил факторизацию колес от колеса 2,3,5, чтобы использовать колесо с фазным коэффициентом 2,3,5,7 с дополнительным предварительным отбором составных чисел, используя дополнительные числа 11, 13 и 17, чтобы значительно уменьшают избыточное отбраковку составных чисел (теперь только отбраковка каждого составного числа в среднем примерно в 1,5 раза). Из-за связанных с DotNet вычислительных накладных расходов (это также относится к колесу 2,3,5 в качестве предыдущей версии) фактическая экономия времени при отбраковке не так уж и велика, но перечисление ответов несколько быстрее из-за многих «простых» составных чисел пропускается в представлении упакованного бита.
  3. Он по-прежнему использует параллельную библиотеку задач (TPL) из DotNet 4 и для многопоточности из пула streamов на основе каждой страницы.
  4. В настоящее время он использует представление базовых простых чисел, которое поддерживает автоматическое расширение массива, содержащегося в этом classе, поскольку в качестве метода streamовой передачи требуется больше базовых простых чисел вместо фиксированного предварительно вычисленного массива базисных чисел, который использовался ранее.
  5. Представление базовых простых чисел было уменьшено до одного байта на базовое число для дальнейшего уменьшения объема памяти; таким образом, общая площадь памяти, отличная от кода, представляет собой массив, который удерживает это представление базовых простых чисел для простых чисел до квадратного корня текущего обрабатываемого диапазона, а буферы с буферами упакованных бит, которые в настоящее время установлены под размером кеша L2 из 256 Килобайтов (наименьший размер страницы 14 586 байт по сравнению с CHNKSZ из 17, как указано) для каждого ядра процессора плюс один дополнительный буфер для обработки переднего плана. С этим кодом около трех мегабайт достаточно для обработки основного диапазона до десяти до четырнадцатой мощности. Помимо скорости, обеспечивающей эффективную многопроцессорную обработку, это уменьшает потребность в памяти, является другим преимуществом использования реализаций постраничного сита.

     class UltimatePrimesSoE : IEnumerable { static readonly uint NUMPRCSPCS = (uint)Environment.ProcessorCount + 1; const uint CHNKSZ = 17; const int L1CACHEPOW = 14, L1CACHESZ = (1 << L1CACHEPOW), MXPGSZ = L1CACHESZ / 2; //for buffer ushort[] //the 2,3,57 factorial wheel increment pattern, (sum) 48 elements long, starting at prime 19 position static readonly byte[] WHLPTRN = { 2,3,1,3,2,1,2,3,3,1,3,2,1,3,2,3,4,2,1,2,1,2,4,3, 2,3,1,2,3,1,3,3,2,1,2,3,1,3,2,1,2,1,5,1,5,1,2,1 }; const uint FSTCP = 11; static readonly byte[] WHLPOS; static readonly byte[] WHLNDX; //to look up wheel indices from position index static readonly byte[] WHLRNDUP; //to look up wheel rounded up index positon values, allow for overfolw static readonly uint WCRC = WHLPTRN.Aggregate(0u, (acc, n) => acc + n); static readonly uint WHTS = (uint)WHLPTRN.Length; static readonly uint WPC = WHTS >> 4; static readonly byte[] BWHLPRMS = { 2, 3, 5, 7, 11, 13, 17 }; const uint FSTBP = 19; static readonly uint BWHLWRDS = BWHLPRMS.Aggregate(1u, (acc, p) => acc * p) / 2 / WCRC * WHTS / 16; static readonly uint PGSZ = MXPGSZ / BWHLWRDS * BWHLWRDS; static readonly uint PGRNG = PGSZ * 16 / WHTS * WCRC; static readonly uint BFSZ = CHNKSZ * PGSZ, BFRNG = CHNKSZ * PGRNG; //number of uints even number of caches in chunk static readonly ushort[] MCPY; //a Master Copy page used to hold the lower base primes preculled version of the page struct Wst { public ushort msk; public byte mlt; public byte xtr; public ushort nxt; } static readonly byte[] PRLUT; /*Wheel Index Look Up Table */ static readonly Wst[] WSLUT; //Wheel State Look Up Table static readonly byte[] CLUT; // a Counting Look Up Table for very fast counting of primes static int count(uint bitlim, ushort[] buf) { //very fast counting if (bitlim < BFRNG) { var addr = (bitlim - 1) / WCRC; var bit = WHLNDX[bitlim - addr * WCRC] - 1; addr *= WPC; for (var i = 0; i < 3; ++i) buf[addr++] |= (ushort)((unchecked((ulong)-2) << bit) >> (i << 4)); } var acc = 0; for (uint i = 0, w = 0; i < bitlim; i += WCRC) acc += CLUT[buf[w++]] + CLUT[buf[w++]] + CLUT[buf[w++]]; return acc; } static void cull(ulong lwi, ushort[] b) { ulong nlwi = lwi; for (var i = 0u; i < b.Length; nlwi += PGRNG, i += PGSZ) MCPY.CopyTo(b, i); //copy preculled lower base primes. for (uint i = 0, pd = 0; ; ++i) { pd += (uint)baseprms[i] >> 6; var wi = baseprms[i] & 0x3Fu; var wp = (uint)WHLPOS[wi]; var p = pd * WCRC + PRLUT[wi]; var pp = (p - FSTBP) >> 1; var k = (ulong)p * (pp + ((FSTBP - 1) >> 1)) + pp; if (k >= nlwi) break; if (k < lwi) { k = (lwi - k) % (WCRC * p); if (k != 0) { var nwp = wp + (uint)((k + p - 1) / p); k = (WHLRNDUP[nwp] - wp) * p - k; if (nwp >= WCRC) wp = 0; else wp = nwp; } } else k -= lwi; var kd = k / WCRC; var kn = WHLNDX[k - kd * WCRC]; for (uint wrd = (uint)kd * WPC + (uint)(kn >> 4), ndx = wi * WHTS + kn; wrd < b.Length; ) { var st = WSLUT[ndx]; b[wrd] |= st.msk; wrd += st.mlt * pd + st.xtr; ndx = st.nxt; } } } static Task cullbf(ulong lwi, ushort[] b, Action f) { return Task.Factory.StartNew(() => { cull(lwi, b); f(b); }); } class Bpa { //very efficient auto-resizing thread-safe read-only indexer class to hold the base primes array byte[] sa = new byte[0]; uint lwi = 0, lpd = 0; object lck = new object(); public uint this[uint i] { get { if (i >= this.sa.Length) lock (this.lck) { var lngth = this.sa.Length; while (i >= lngth) { var bf = (ushort[])MCPY.Clone(); if (lngth == 0) { for (uint bi = 0, wi = 0, w = 0, msk = 0x8000, v = 0; w < bf.Length; bi += WHLPTRN[wi++], wi = (wi >= WHTS) ? 0 : wi) { if (msk >= 0x8000) { msk = 1; v = bf[w++]; } else msk <<= 1; if ((v & msk) == 0) { var p = FSTBP + (bi + bi); var k = (p * p - FSTBP) >> 1; if (k >= PGRNG) break; var pd = p / WCRC; var kd = k / WCRC; var kn = WHLNDX[k - kd * WCRC]; for (uint wrd = kd * WPC + (uint)(kn >> 4), ndx = wi * WHTS + kn; wrd < bf.Length; ) { var st = WSLUT[ndx]; bf[wrd] |= st.msk; wrd += st.mlt * pd + st.xtr; ndx = st.nxt; } } } } else { this.lwi += PGRNG; cull(this.lwi, bf); } var c = count(PGRNG, bf); var na = new byte[lngth + c]; sa.CopyTo(na, 0); for (uint p = FSTBP + (this.lwi << 1), wi = 0, w = 0, msk = 0x8000, v = 0; lngth < na.Length; p += (uint)(WHLPTRN[wi++] << 1), wi = (wi >= WHTS) ? 0 : wi) { if (msk >= 0x8000) { msk = 1; v = bf[w++]; } else msk <<= 1; if ((v & msk) == 0) { var pd = p / WCRC; na[lngth++] = (byte)(((pd - this.lpd) << 6) + wi); this.lpd = pd; } } this.sa = na; } } return this.sa[i]; } } } static readonly Bpa baseprms = new Bpa(); static UltimatePrimesSoE() { WHLPOS = new byte[WHLPTRN.Length + 1]; //to look up wheel position index from wheel index for (byte i = 0, acc = 0; i < WHLPTRN.Length; ++i) { acc += WHLPTRN[i]; WHLPOS[i + 1] = acc; } WHLNDX = new byte[WCRC + 1]; for (byte i = 1; i < WHLPOS.Length; ++i) { for (byte j = (byte)(WHLPOS[i - 1] + 1); j <= WHLPOS[i]; ++j) WHLNDX[j] = i; } WHLRNDUP = new byte[WCRC * 2]; for (byte i = 1; i < WHLRNDUP.Length; ++i) { if (i > WCRC) WHLRNDUP[i] = (byte)(WCRC + WHLPOS[WHLNDX[i - WCRC]]); else WHLRNDUP[i] = WHLPOS[WHLNDX[i]]; } Func nmbts = (v) => { var acc = 0; while (v != 0) { acc += (int)v & 1; v >>= 1; } return acc; }; CLUT = new byte[1 << 16]; for (var i = 0; i < CLUT.Length; ++i) CLUT[i] = (byte)nmbts((ushort)(i ^ -1)); PRLUT = new byte[WHTS]; for (var i = 0; i < PRLUT.Length; ++i) { var t = (uint)(WHLPOS[i] * 2) + FSTBP; if (t >= WCRC) t -= WCRC; if (t >= WCRC) t -= WCRC; PRLUT[i] = (byte)t; } WSLUT = new Wst[WHTS * WHTS]; for (var x = 0u; x < WHTS; ++x) { var p = FSTBP + 2u * WHLPOS[x]; var pr = p % WCRC; for (uint y = 0, pos = (p * p - FSTBP) / 2; y < WHTS; ++y) { var m = WHLPTRN[(x + y) % WHTS]; pos %= WCRC; var posn = WHLNDX[pos]; pos += m * pr; var nposd = pos / WCRC; var nposn = WHLNDX[pos - nposd * WCRC]; WSLUT[x * WHTS + posn] = new Wst { msk = (ushort)(1 << (int)(posn & 0xF)), mlt = (byte)(m * WPC), xtr = (byte)(WPC * nposd + (nposn >> 4) - (posn >> 4)), nxt = (ushort)(WHTS * x + nposn) }; } } MCPY = new ushort[PGSZ]; foreach (var lp in BWHLPRMS.SkipWhile(p => p < FSTCP)) { var p = (uint)lp; var k = (p * p - FSTBP) >> 1; var pd = p / WCRC; var kd = k / WCRC; var kn = WHLNDX[k - kd * WCRC]; for (uint w = kd * WPC + (uint)(kn >> 4), ndx = WHLNDX[(2 * WCRC + p - FSTBP) / 2] * WHTS + kn; w < MCPY.Length; ) { var st = WSLUT[ndx]; MCPY[w] |= st.msk; w += st.mlt * pd + st.xtr; ndx = st.nxt; } } } struct PrcsSpc { public Task tsk; public ushort[] buf; } class nmrtr : IEnumerator, IEnumerator, IDisposable { PrcsSpc[] ps = new PrcsSpc[NUMPRCSPCS]; ushort[] buf; public nmrtr() { for (var s = 0u; s < NUMPRCSPCS; ++s) ps[s] = new PrcsSpc { buf = new ushort[BFSZ] }; for (var s = 1u; s < NUMPRCSPCS; ++s) { ps[s].tsk = cullbf((s - 1u) * BFRNG, ps[s].buf, (bfr) => { }); } buf = ps[0].buf; } ulong _curr, i = (ulong)-WHLPTRN[WHTS - 1]; int b = -BWHLPRMS.Length - 1; uint wi = WHTS - 1; ushort v, msk = 0; public ulong Current { get { return this._curr; } } object IEnumerator.Current { get { return this._curr; } } public bool MoveNext() { if (b < 0) { if (b == -1) b += buf.Length; //no yield!!! so automatically comes around again else { this._curr = (ulong)BWHLPRMS[BWHLPRMS.Length + (++b)]; return true; } } do { i += WHLPTRN[wi++]; if (wi >= WHTS) wi = 0; if ((this.msk <<= 1) == 0) { if (++b >= BFSZ) { b = 0; for (var prc = 0; prc < NUMPRCSPCS - 1; ++prc) ps[prc] = ps[prc + 1]; ps[NUMPRCSPCS - 1u].buf = buf; ps[NUMPRCSPCS - 1u].tsk = cullbf(i + (NUMPRCSPCS - 1u) * BFRNG, buf, (bfr) => { }); ps[0].tsk.Wait(); buf = ps[0].buf; } v = buf[b]; this.msk = 1; } } while ((v & msk) != 0u); _curr = FSTBP + i + i; return true; } public void Reset() { throw new Exception("Primes enumeration reset not implemented!!!"); } public void Dispose() { } } public IEnumerator GetEnumerator() { return new nmrtr(); } IEnumerator IEnumerable.GetEnumerator() { return new nmrtr(); } static void IterateTo(ulong top_number, Action actn) { PrcsSpc[] ps = new PrcsSpc[NUMPRCSPCS]; for (var s = 0u; s < NUMPRCSPCS; ++s) ps[s] = new PrcsSpc { buf = new ushort[BFSZ], tsk = Task.Factory.StartNew(() => { }) }; var topndx = (top_number - FSTBP) >> 1; for (ulong ndx = 0; ndx <= topndx; ) { ps[0].tsk.Wait(); var buf = ps[0].buf; for (var s = 0u; s < NUMPRCSPCS - 1; ++s) ps[s] = ps[s + 1]; var lowi = ndx; var nxtndx = ndx + BFRNG; var lim = topndx < nxtndx ? (uint)(topndx - ndx + 1) : BFRNG; ps[NUMPRCSPCS - 1] = new PrcsSpc { buf = buf, tsk = cullbf(ndx, buf, (b) => actn(lowi, lim, b)) }; ndx = nxtndx; } for (var s = 0u; s < NUMPRCSPCS; ++s) ps[s].tsk.Wait(); } public static long CountTo(ulong top_number) { if (top_number < FSTBP) return BWHLPRMS.TakeWhile(p => p <= top_number).Count(); var cnt = (long)BWHLPRMS.Length; IterateTo(top_number, (lowi, lim, b) => { Interlocked.Add(ref cnt, count(lim, b)); }); return cnt; } public static ulong SumTo(uint top_number) { if (top_number < FSTBP) return (ulong)BWHLPRMS.TakeWhile(p => p <= top_number).Aggregate(0u, (acc, p) => acc += p); var sum = (long)BWHLPRMS.Aggregate(0u, (acc, p) => acc += p); Func sumbf = (lowi, bitlim, buf) => { var acc = 0L; for (uint i = 0, wi = 0, msk = 0x8000, w = 0, v = 0; i < bitlim; i += WHLPTRN[wi++], wi = wi >= WHTS ? 0 : wi) { if (msk >= 0x8000) { msk = 1; v = buf[w++]; } else msk <<= 1; if ((v & msk) == 0) acc += (long)(FSTBP + ((lowi + i) << 1)); } return acc; }; IterateTo(top_number, (pos, lim, b) => { Interlocked.Add(ref sum, sumbf(pos, lim, b)); }); return (ulong)sum; } static void IterateUntil(Func prdct) { PrcsSpc[] ps = new PrcsSpc[NUMPRCSPCS]; for (var s = 0u; s < NUMPRCSPCS; ++s) { var buf = new ushort[BFSZ]; ps[s] = new PrcsSpc { buf = buf, tsk = cullbf(s * BFRNG, buf, (bfr) => { }) }; } for (var ndx = 0UL; ; ndx += BFRNG) { ps[0].tsk.Wait(); var buf = ps[0].buf; var lowi = ndx; if (prdct(lowi, buf)) break; for (var s = 0u; s < NUMPRCSPCS - 1; ++s) ps[s] = ps[s + 1]; ps[NUMPRCSPCS - 1] = new PrcsSpc { buf = buf, tsk = cullbf(ndx + NUMPRCSPCS * BFRNG, buf, (bfr) => { }) }; } } public static ulong ElementAt(long n) { if (n < BWHLPRMS.Length) return (ulong)BWHLPRMS.ElementAt((int)n); long cnt = BWHLPRMS.Length; var ndx = 0UL; var cycl = 0u; var bit = 0u; IterateUntil((lwi, bfr) => { var c = count(BFRNG, bfr); if ((cnt += c) < n) return false; ndx = lwi; cnt -= c; c = 0; do { var w = cycl++ * WPC; c = CLUT[bfr[w++]] + CLUT[bfr[w++]] + CLUT[bfr[w]]; cnt += c; } while (cnt < n); cnt -= c; var y = (--cycl) * WPC; ulong v = ((ulong)bfr[y + 2] << 32) + ((ulong)bfr[y + 1] << 16) + bfr[y]; do { if ((v & (1UL << ((int)bit++))) == 0) ++cnt; } while (cnt <= n); --bit; return true; }); return FSTBP + ((ndx + cycl * WCRC + WHLPOS[bit]) << 1); } } 

The above code takes about 59 milliseconds to find the primes to two million (slightly slower than some of the other simpler codes due to initialization overhead), but calculates the primes to one billion and the full number range in 1.55 and 5.95 seconds, respectively. This isn't much faster than the last version due to the DotNet extra overhead of an extra array bound check in the enumeration of found primes compared to the time expended in culling composite numbers which is less than a third of the time spent emumerating, so the saving in culling composites is cancelled out by the extra time (due to an extra array bounds check per prime candidate) in the enumeration. However, for many tasks involving primes, one does not need to enumerate all primes but can just compute the answers without enumeration.

For the above reasons, this class provides the example static methods "CountTo", "SumTo", and "ElementAt" to count or sum the primes to a given upper limit or to output the zero-based nth prime, respectively. The "CountTo" method will produce the number of primes to one billion and in the 32-bit number range in about 0.32 and 1.29 seconds, respectively ; the "ElementAt" method will produce the last element in these ranges in about 0.32 and 1.25 seconds, respectively , and the "SumTo" method produces the sum of all the primes in these ranges in about 0.49 and 1.98 seconds respectively . This program calculates the sum of all the prime numbers to four billion plus as here in less time than many naive implementations can sum all the prime numbers to two million as in Euler Problem 10, for over 2000 times the practical range!

This code is only about four times slower than very highly optimized C code used by primesieve , and the reasons it is slower are mostly due to DotNet, as follows (discussing the case of a 256 Kilobyte buffer, which is the size of the L2 cache):

  1. Most of the execution time is spent in the main composite culling loop, which is the last "for loop" in the private static "cull" method" and only contains four statements per loop plus the range check.
  2. In DotNet, this compiles to take about 21.83 CPU clock cycles per loop, including about 5 clock cycles for the two array bounds checks per loop.
  3. The very efficient C compiler converts this loop into only about 8.33 clock cycles for an advantage of about 2.67 times.
  4. Primesieve also uses extreme manual "loop unrolling" optimizations to reduce the average time to perform the work per loop to about 4.17 clock cycles per composite cull loop, for a additional gain of two times and a total gain of about 5.3 times.
  5. Now, the highly optimized C code both doesn't Hyper Thread (HT) as well as the less efficient Just In Time (JIT) compiler produced code and as well, the OpemMP multi-threading used by primesieve doesn't appear to be as well adapted to this problem as use of the Thread Pool threads here, so the final multi-threaded gain in only about four times.
  6. One might consider the use of "unsafe" pointers to eliminate the array bounds checks, and it was tried, but the JIT compiler does not optimize pointers as well as normal array based code, so the gain of not having array bound checks is cancelled by the less efficient code (every pointer access (re)loads the pointer address from memory rather than using a register already pointing to that address as in the optimized array case).
  7. Primesieve is even faster when using smaller buffer sizes as in the size of the available L1 cache (16 Kilobytes when multi-threading for the i3/i5/i7 CPU's) as its more efficient code has more of an advantage reducing the average memory access time to one clock cycle from about four clock cycles, which advantage makes much less of a difference to the DotNet code, which gains more from less processing per reduced number of pages. Thus, primesieve is about five times faster when each use their most efficient buffer size.

This DotNet code will count (CountTo) the number of primes to ten to the power thirteen (ten trillion) in about an hour and a half (tested) and the number of primes to a hundred trillion (ten to the fourteenth) in just over a half day (estimated) as compared to twenty minutes and under four hours for primesieve, respectively. This is interesting historically as until 1985 only the count of primes in the range to ten to the thirteenth was known since it would take too much execution time on the (expensive) supercomputers of that day to find the range ten times larger; now we can easily compute the number of primes in those ranges on a common desktop computer (in this case, an Intel i7-2700K - 3.5 GHz)!

Using this code, it is easy to understand why Professor Atkin and Bernstein thought that the SoA was faster than the SoE - a myth that persists to this day , with the reasoning as follows:

  1. It is easy to have any SoA implementation count the number of state toggles and square free composite number culls (which last can be optimized using the same 2,3,5 wheel optimization as used inherently by the SoA algorithm) to determine that the total number of both of these operations is about 1.4 billion for the 32-bit number range.
  2. Bernstein's "equivalent" SoE implementation to his SoA implementation (neither of which are as optimized as this code), which uses the same 2,3,5 wheel optimization as the SoA, will have a total of about 1.82 billion cull operations with the same cull loop computational complexity.
  3. Therefore, Bernstein's results of about 30% better performance as compared to his implementation of the SoE is about right, just based on the number of equivalent operations. However, his implementation of the SoE did not take wheel factorization "to the max" since the SoA does not much respond to further degrees of wheel factorization as the 2,3,5 wheel is "baked in" to the basic algorithm.
  4. The wheel factorization optimizations used here reduces the number of composite cull operations to about 1.2 billion for the 32-bit number range; therefore, this algorithm using this degree of wheel factorization will run about 16.7% faster than an equivalent version of SoA since the culling loop(s) can be implemented about the same for each algorithm.
  5. The SoE with this level of optimizations is easier to write than an equivalent SoA as there needs only be one state table look up array to cull across the base primes instead of an additional state look up array's for each of the four quadratic equation solution cases that produce valid state toggles.
  6. When written using equivalent implementations as this code, the code will respond equivalently to C compiler optimizations for SoA as for the SoE used in primesieve.
  7. The SoA implementation will also respond to the extreme manual "loop unrolling" optimization used by primesieve just as effectively as for the SoE implementation.
  8. Therefore, if I went though all the work to implement the SoA algorithm using the techniques as for the above SoE code, the resulting SoA would only be a slight amount slower when the output primes were enumerated but would be about 16.7% slower when using the static direct methods .
  9. Memory footprint of the two algorithms is no different as both require the representation of the base primes and the same number of segment page buffers.

EDIT_ADD: Interestingly, this code runs 30% faster in x86 32-bit mode than in x64 64-bit mode, likely due to avoiding the slight extra overhead of extending the uint32 numbers to ulong's. All of the above timings are for 64-bit mode. END_EDIT_ADD

In summary: The segment paged Sieve of Atkin is actually slower than a maximally optimized segment paged Sieve of Eratosthenes with no saving in memory requirements!!!

I say again: "Why use the Sieve of Atkin?" ,

This is a more explicit answer to the question(s) raised, as follows:

Is there a way (the incremental Sieve of Atkin – GBG) can be done?

Of course the Sieve of Atkin (SoA) can be implemented as a segmented algorithm and in fact does not require the use of segmenting arrays as suggested at all but rather can be done just using raw sequences as I have done functionally using F# , although this is quite a bit slower than the extra efficiency permitted by using mutable arrays for the culling.

I am thinking it could go something like this: 1. Start with some trivial limit 2. Find all the primes up to the limit 3. Find all the primes up to the limit 4. Yield all the newfound primes 5. Increase the limit (by doubling or squaring the old limit or something like that) 6. Goto step 2

The above algorithm can be implemented in at least two different ways: 1) save the state of ‘x’ and ‘y’ for each value of ‘x’ when the sequences “run off” the current segment and start up again with those values for the next segment, or 2) calculate the applicable pair values of ‘x’ and ‘y’ to use for the new segment. Although the first way is simpler, I recommend the second method for two reasons: 1) it doesn’t use memory for all the (many) values of x and y that must be saved and only a representation of the base primes must be saved in memory for the “squares free” culling step, and 2) it opens the door to using multi-threading and assigning independent thread operations for each page segment for large savings in time on a multi-processor computer.

And indeed, a better understanding of ‘x’ and ‘y’ is necessary:

My main problem is that I don’t quite understand what x and y for example is in this algorithm. Like, could I just use the same algorithm kind of but set x and y to oldLimit (initially 1) and run it up to newLimit?

There has been one answer addressing this but perhaps it isn’t explicit enough. It may be easier to think of these quadratic equations as a potentially infinite sequence of sequences where one of ‘x’ or ‘y’ is fixed starting from their lowest values and the other variable produces all of the elements per sequence. For instance, one could consider the odd expression of the quadratic equation “4*x^2 + y^2” as the sequence of sequences starting with 5, 17, 37, 65, …, and with each of these sequences have elements as in {5, 13, 29, 53, …}, {17, 25, 41, 65, …}, {37, 45, 61, 85, …}, {65, 73, 89, 125, …}, … Obviously, some of these elements are not prime as they are composites of 3 or 5 and that is why these need to be eliminated by either the modulo test, or alternatively as in the Bernstein implementation, these can be skipped automatically by recognizing patterns in the modulo arithmetic for the generators so they never actually even appear.

Implementing the first easier way of producing a segmented version of the SoA then just requires saving the state of each of the sequence of sequences, which is basically what is done in the F# incremental implementation (although using a folded tree structure for efficiency), which could easily be adapted to working over an array page range. For the case where the sequence state is computed at the beginning of every segment page, one just needs to compute how many elements would fit into the space up to the number represented by the lowest element in the new segment page for each “active” sequence, where “active” means a sequence whose starting element is lower than the number represented by the start index of the segment page.

As for pseudo-code on how to implement the segmentation of an array based SoA sieve, I have written something up for a related post , which shows how this can be accomplished.

The point of this is so that I won’t have to set that limit. So that I can for example use Linq and just Take() however many primes I need, not worrying about if the limit is high enough, et cetera.

As stated in another answer, you could accomplish this end goal by just setting the maximum “limit” as a constant in your code, but this would be quite inefficient for small ranges of primes as culling would be taking place over a much larger array than necessary. As already stated, other than for better efficiency and reducing memory use by a huge factor, segmentation also has other advantages in permitting the efficient use of multi-processing. However, use of the Take(), TakeWhile(), Where(), Count(), etcetera, methods won’t provide very good performance for large ranges of primes as their use involves many stacked method calls per element at many clock cycles per call and return. But you would have the option of either using these or more imperative forms of program flow, so this isn’t a real objection.

I can try to explain what x and y does, but I don’t think you can do what you ask without restarting the loops from the beginning. This is pretty much the same for any “sieve”-algorithm.

What the sieve does is basically count how many different quadratic equations (even or odd) have each number as a solution. The specific equation checked for each number is different depending on what n % 12 is.

For example, numbers n which have a mod 12 remainder of 1 or 5 are prime if and only if the number of solutions to 4* x ^2+ y ^2= n is odd and the number is square-free. The first loop simply loops through all possible values of x and y that could satisfy these different equations. By flipping the isPrime[ n ] each time we find a solution for that n , we can keep track of whether the number of solutions is odd or even.

The thing is that we count this for all possible n at the same time, which makes this much more efficient than checking one at the time. Doing it for just some n would take more time (because you would have to make sure that n >= lower_limit in the first loop) and complicate the second loop, since that one requires knowledge of all primes smaller than sqrt.

The second loop checks that the number is square-free (has no factor which is a square of a prime number).

Also, I don’t think your implementation of the sieve of Atkin is necessarily faster than a straight-forward sieve of Eratosthenes would be. However, the fastest possible most optimized sieve of Atkin would beat the fastest possible most optimized sieve of Eratosthenes.

Давайте будем гением компьютера.