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

Ответ 1

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

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

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

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

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

static IEnumerable<uint> 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 MegaBytes оперативной памяти).

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

static IEnumerable<uint> 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<uint> 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 процессора, ускоряются из-за более эффективного доступа к памяти и во-вторых, поскольку сегментация позволяет легко разбивать задание на страницы, которые можно отбирать в фоновом режиме, используя многоядерные процессоры для ускорения, которые могут быть пропорциональны количеству используемых сердечников.

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

static IEnumerable<uint> 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<BUFSZ; k+=pa[WHLPTRN[sw]-1], sw=(sw<7u) ? ++sw : 0u) buf[(int)k]=false; } }
    if (buf[(int)si]) yield return 7u + i + i; } }

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

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

EDIT_ADD:

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

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

class fastprimesSoE : IEnumerable<uint>, 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<uint>, 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<uint> GetEnumerator() { return new nmrtr(); }
  IEnumerator IEnumerable.GetEnumerator() { return new nmrtr(); } }

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

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

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

Ответ 2

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

• Диапазон использования был увеличен до 64-но беззнакового числа диапазон 18,446,744,073,709,551,615 с проверками переполнения диапазона удалено, так как маловероятно, что вы захотите запустить программу в течение сотен лет, которые потребовались бы для обработки всего спектра чисел до этого предела. Это связано с очень небольшими затратами на обработку так как пейджинг может быть выполнен с использованием 32-битных диапазонов страниц и только конечный первичный вывод должен быть вычислен как 64-разрядное число.
• Он увеличил факторизацию колес от колеса 2,3,5, чтобы использовать 2,3,5,7 колесо фаз-фактора с дополнительным предварительным отбором композитных числа с использованием дополнительных простых чисел 11, 13 и 17, чтобы значительно уменьшить избыточное отбракование составных чисел (теперь только отбраковывая каждое составное число в среднем примерно в 1,5 раза). В связи (связанные с DotNet) вычислительные накладные расходы на это (также относится к колесу 2,3,5 в качестве предыдущей версии) фактическое время сохранение в отбраковке не так уж и велико, но перечисление ответов несколько быстрее из-за многих "простых" составных чисел пропущен в представлении упакованного бита.
• Он по-прежнему использует параллельную библиотеку задач (TPL) от DotNet 4 и выше для многопоточности из пула потоков на основе каждой страницы.
• Теперь он использует представление базовых простых чисел, которое автоматически поддерживает увеличивая массив, содержащийся в этом классе, поскольку больше базовых простых чисел требуемый как метод безопасности потока, а не фиксированный предварительно вычисленный массив базовых чисел, используемый ранее.

• Представление базовых простых чисел сокращено до одного байта на базу prime для дальнейшего уменьшения объема памяти; таким образом, общая сумма область памяти, отличная от кода, - это массив для хранения этой базы Представление простых чисел для простых чисел до квадратного корня из обрабатываемый текущий диапазон, и буферы буферов упакованных бит, которые в настоящее время установлены под размером кеша L2, равным 256 килобайтам (наименьший размер страницы составляет 14 586 байтов, а CHNKSZ - 17 ) для каждого ядра процессора плюс один дополнительный буфер для переднего плана задача для обработки. С этим кодом около трех мегабайт достаточный для обработки основного диапазона до десяти-четырнадцатого мощность. Помимо скорости благодаря эффективной многопроцессорной обработке, это сокращение потребности в памяти является другим преимуществом использования paging sieve.

  class UltimatePrimesSoE : IEnumerable<ulong> {
    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<ushort[]> 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<ushort, int> 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<ulong>, 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<ulong> GetEnumerator() { return new nmrtr(); }
    IEnumerator IEnumerable.GetEnumerator() { return new nmrtr(); }
    static void IterateTo(ulong top_number, Action<ulong, uint, ushort[]> 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<ulong, uint, ushort[], long> 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<ulong, ushort[], bool> 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); } }
  

Приведенный выше код занимает около 59 миллисекунд, чтобы найти простые числа до двух миллионов (немного медленнее, чем некоторые другие простые коды из-за издержек инициализации), но вычисляет простые числа до одного миллиарда и полный диапазон номеров в 1,55 и 5,95 секунды, соответственно. Это не намного быстрее, чем последняя версия из-за дополнительных накладных расходов DotNet дополнительной проверки привязки массива в перечислении найденных простых чисел по сравнению с временем, затраченным на отбор составных чисел, составляющих менее трети времени, потраченного на эмумеризацию, поэтому сохранение в отбраковке композитов отменяется дополнительным временем (из-за дополнительной проверки границ массива на первого кандидата) в перечислении. Однако для многих задач с участием простых чисел не нужно перечислять все простые числа, но может просто вычислять ответы без перечисления.

По вышеуказанным причинам этот класс предоставляет пример статических методов "CountTo", "SumTo" и "ElementAt" для подсчета или суммирования простых чисел до заданного верхнего предела или для вывода нулевого числа n-го числа соответственно. Метод "CountTo" приведет к количеству простых чисел до одного миллиарда и в 32-битном диапазоне номеров примерно в 0,32 и 1,29 секунды, соответственно > ; метод "ElementAt" будет производить последний элемент в этих диапазонах примерно 0,32 и 1,25 секунды, соответственно, а метод "SumTo сумма всех простых чисел в этих диапазонах примерно 0,49 и 1,98 секунды соответственно. Эта программа вычисляет сумму всех простых чисел до четырех миллиардов плюс, как здесь за меньшее время, чем многие наивные реализации могут суммировать все простые числа до двух миллионов, как в Эйлере 10, более чем в 2000 раз больше практического диапазона!

Этот код работает примерно в четыре раза медленнее, чем очень высоко оптимизированный код C, используемый primesieve, а причины этого медленнее, в основном из-за DotNet, следующим образом (обсуждение вопроса о 256-килобайтном буфере, который является размером кэша L2):

• Большая часть времени выполнения расходуется на основной сборке loop, который является последним "для цикла" в частной статической "cull", метод "и содержит только четыре оператора для каждого цикла плюс диапазон проверить.
• В DotNet этот компилятор принимает около 21,83 тактовых такта ЦП на каждый цикл, включающий около 5 тактов для двух границ массива проверяет на цикл.
• Очень эффективный компилятор C преобразует этот цикл только в 8,33 тактовых цикла для преимущества около 2,67 раза.
• Primesieve также использует экстремальные ручные "разворачивание петли" для сократить среднее время выполнения работы за цикл до примерно 4,17 тактовых циклов на композитный цикл, для дополнительного усиления двух раз и общий прирост примерно в 5,3 раза.
• Теперь высоко оптимизированный код C не является Hyper Thread (HT) как так как менее эффективный компилятор Just In Time (JIT) выпустил код и, кроме того, многопоточность OpemMP, используемая primesieve, не как представляется, также адаптированы к этой проблеме, поскольку использование Thread Объемы потоков здесь, поэтому окончательный многопоточный выигрыш в четыре раза.
• Можно подумать об использовании "небезопасных" указателей для устранения проверка границ массива, и она была проверена, но компилятор JIT не оптимизировать указатели, а также обычный код на основе массива, поэтому коэффициент усиления отсутствие проверки границ массива отменяется менее эффективным код (каждый доступ указателя (re) загружает адрес указателя из памяти вместо использования регистра, уже указывающего на этот адрес, как в оптимизированный массив).
• Primesieve работает быстрее при использовании меньших размеров буфера, как в размер доступного кеша L1 (16 килобайт при многопоточности для процессоров i3/i5/i7), поскольку его более эффективный код имеет больше преимущество уменьшает среднее время доступа к памяти на один такт от примерно четырех тактовых циклов, преимущество которых значительно меньше разница с кодом DotNet, который получает больше от меньшего обработки за меньшее количество страниц. Таким образом, в пять раз быстрее, когда каждый использует их наиболее эффективный размер буфера.

Этот код DotNet будет подсчитывать (CountTo) количество простых чисел до десяти до тринадцати (десять триллионов) в течение примерно полутора часов (проверенных) и количество простых чисел до ста триллионов (от десяти до четырнадцатого) в течение чуть более полудня (по оценкам), по сравнению с двадцатью минутами и менее чем за четыре часа соответственно. Это исторически интересно, так как до 1985 года было известно только количество простых чисел в диапазоне от десяти до тринадцатого, поскольку для суперкомпьютеров этого дня потребовалось бы слишком много времени, чтобы найти диапазон в десять раз больше; теперь мы можем легко вычислить количество простых чисел в этих диапазонах на общем настольном компьютере (в данном случае Intel i7-2700K - 3,5 ГГц)!

Используя этот код, легко понять, почему профессор Аткин и Бернштейн считали, что SoA быстрее, чем SoE - миф, который сохраняется и по сей день, с рассуждением следующим образом:

• Легко получить, что любая реализация SoA подсчитывает количество состояний переключатели и квадратные квадратные отборные числовые отрезки (последние могут быть оптимизирована с использованием той же оптимизации 2,3,5 колес, которая используется по существу по алгоритму SoA), чтобы определить, что общее число обоих эти операции составляют около 1,4 млрд. для 32-битного диапазона номеров.
• Бернштейн "эквивалент" реализации SoE для его SoA (ни одна из них не оптимизирована, как этот код), который использует ту же 2,3,5-колесную оптимизацию, что и SoA, будет иметь всего около 1,82 млрд. операций cull с одной и той же схемой отрезка вычислительная сложность.
• Таким образом, результаты Бернштейна примерно на 30% лучше, чем по сравнению с его реализация SoE - это правильно, просто основанный на количестве эквивалентных операций. Однако его внедрение SoE не привело к факторизации колес "max", так как SoA не очень реагирует на дополнительные градусы колеса факторизация, так как колесо 2,3,5 "запекается" в основной алгоритм.
• Используемые здесь оптимизации факторизации колес уменьшают количество составные операции отбирания примерно до 1,2 млрд для 32-разрядного номера ассортимент; поэтому этот алгоритм с использованием этой степени колеса факторизация будет работать на 16,7% быстрее, чем эквивалентная версия от SoA, так как цикл отсечения может быть реализован примерно одинаково для каждый алгоритм.
• SoE с таким уровнем оптимизации легче писать, чем эквивалент SoA, так как нужен только один массив поиска таблицы состояний отбирать базовые простые числа вместо дополнительного внешнего вида для каждого из четырех решений уравнения квадратичного уравнения которые производят правильные переключения состояний.
• При написании с использованием эквивалентных реализаций в качестве этого кода код будет отвечать эквивалентно оптимизациям компилятора C для SoA, так как для SoE используется в primesieve.
• Реализация SoA также будет реагировать на экстремальный ручной цикл "разворачивание", оптимизированное так же эффективно, как и для реализации SoE.
• Поэтому, если бы я пошел, хотя вся работа по внедрению SoA алгоритм с использованием методов, как для вышеуказанного кода SoE, в результате SoA будет только немного меньше, когда выход простые числа были перечислены, но будут примерно на 16,7% медленнее при использовании статические прямые методы.
• Объем памяти двух алгоритмов не отличается, так как оба требуют представления базовых простых чисел и того же числа сегментных буферов страницы.

EDIT_ADD: Интересно, что этот код работает на 30% быстрее в 32-разрядном режиме x86, чем в 64-битном режиме 64, вероятно, из-за того, что он избегает незначительных дополнительных накладных расходов, расширяя номера uint32 до ulong. Все приведенные выше тайминги предназначены для 64-битного режима. END_EDIT_ADD

В заключение: сегмент paged Сито Аткина на самом деле медленнее, чем максимально оптимизированный сегмент, выгружаемый ситом Eratosthenes без необходимости сохранения в памяти!

Я снова говорю: "Зачем использовать сито Аткина?" .

Ответ 3

Это более явный ответ на вопрос (ы), поднятый следующим образом:

Есть ли способ (инкрементное сито Аткина - GBG)?

Конечно, сито Atkin (SoA) может быть реализовано как сегментированный алгоритм и на самом деле не требует использования сегментированных массивов, как предлагается вообще, а может быть сделано только с использованием необработанных последовательностей как Я сделал функционально, используя F #, хотя это довольно немного медленнее, чем дополнительная эффективность, разрешенная с помощью изменяемых массивов для отбраковки.

Я думаю, что это может пойти примерно так:  1. Начните с некоторого тривиального предела  2. Найти все простые числа до предела  3. Найдите все простые числа до предела  4. Допустим все новообретенные простые числа  5. Увеличьте предел (удвоив или возведя в квадрат старый лимит или что-то в этом роде)  6. Перейдите к шагу 2

Вышеупомянутый алгоритм может быть реализован, по крайней мере, двумя разными способами: 1) сохранить состояние "x" и "y" для каждого значения "x", когда последовательности "выбегут" из текущего сегмента и снова запустите с этими значениями для следующего сегмента или 2) рассчитать применимые значения пары 'x' и 'y' для использования для нового сегмента. Хотя первый способ проще, я рекомендую второй метод по двум причинам: 1) он не использует память для всех (многих) значений x и y, которые должны быть сохранены, и должно быть сохранено только представление базовых простых чисел в памяти для шага отсечения "без квадратов" и 2) он открывает дверь для использования многопоточности и назначения независимых операций потока для каждого сегмента страницы для значительной экономии времени на многопроцессорном компьютере.

И действительно, лучше понимать "х" и "у":

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

Был один ответ, касающийся этого, но, возможно, он недостаточно ясен. Может быть, проще думать об этих квадратичных уравнениях как о потенциально бесконечной последовательности последовательностей, где один из "х" или "у" фиксирован, начиная с их младших значений, а другая переменная производит все элементы на каждую последовательность. Например, можно было бы рассматривать нечетное выражение квадратичного уравнения "4 * x ^ 2 + y ^ 2" как последовательность последовательностей, начинающихся с 5, 17, 37, 65,... и каждая из этих последовательностей имеет элементы, как в {5, 13, 29, 53,...}, {17, 25, 41, 65,...}, {37, 45, 61, 85,...}, {65, 73, 89, 125,...},... Очевидно, что некоторые из этих элементов не являются первичными, так как они представляют собой композиты из 3 или 5, и поэтому их необходимо устранить либо по модулю, либо в качестве альтернативы, как в Bernstein они могут быть пропущены автоматически, распознавая шаблоны в модульной арифметике для генераторов, чтобы они никогда даже не отображались.

Реализация первого более простого способа создания сегментированной версии SoA требует просто сохранения состояния каждой из последовательностей последовательностей, что в основном происходит в инкрементной реализации F # (хотя использование сложенной древовидной структуры для эффективности), который можно легко адаптировать для работы над диапазоном страниц массива. В случае, когда состояние последовательности вычисляется в начале каждой страницы сегмента, нужно просто вычислить, сколько элементов будет вписываться в пространство до числа, представленного самым низким элементом на странице нового сегмента для каждой "активной" последовательности, где "active" означает последовательность, стартовый элемент которой меньше числа, представленного начальным индексом страницы сегмента.

Что касается псевдокода о том, как реализовать сегментацию массива SoA, я написал что-то для связанного сообщения, которое показывает, как это может выполняются.

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

Как указано в другом ответе, вы можете достичь этой конечной цели, просто установив максимальный "предел" в качестве константы в вашем коде, но это было бы довольно неэффективно для небольших диапазонов простых чисел, поскольку отбраковка будет происходить на протяжении многих чем требуется. Как уже было сказано, кроме повышения эффективности и сокращения использования памяти огромным фактором, сегментация также имеет другие преимущества в разрешении эффективного использования многопроцессорной обработки. Однако использование методов Take(), TakeWhile(), Where(), Count() и т.д. Не обеспечит очень хорошую производительность для больших диапазонов простых чисел, так как их использование включает в себя множество вызовов вызовов в стеке для каждого элемента на многих тактах за звонок и возврат. Но у вас будет возможность использовать эти или более императивные формы потока программы, поэтому это не является реальным возражением.

Ответ 4

Я могу попытаться объяснить, что делают x и y, но я не думаю, что вы можете делать то, что вы просите, без перезапуска циклов с самого начала. Это почти то же самое для любого "ситового" алгоритма.

Что такое сито, в основном подсчитывает, сколько разных квадратичных уравнений (четных или нечетных) имеют каждое число в качестве решения. Конкретное уравнение, проверенное для каждого числа, отличается в зависимости от того, что такое n% 12.

Например, числа n, которые имеют остаток mod 12 от 1 или 5, являются первичными тогда и только тогда, когда число решений для 4 * x ^ 2 + y ^ 2 = n нечетно, а число является квадратным. Первый цикл просто пересекает все возможные значения x и y, которые могут удовлетворять этим различным уравнениям. Перебрасывая isPrime [n] каждый раз, когда мы находим решение для этого n, мы можем отслеживать, является ли число решений нечетным или четным.

Дело в том, что мы подсчитываем это для всех возможных n одновременно, что делает его намного более эффективным, чем проверка на тот момент. Выполнение этого только для некоторого n потребует больше времени (потому что вам нужно будет убедиться, что n >= lower_limit в первом цикле) и усложнить второй цикл, поскольку для этого требуется знание всех простых чисел, меньших, чем sqrt.

Второй цикл проверяет, что число является квадратным (не имеет коэффициента, который является квадратом простого числа).

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