Алгоритм нахождения N первых простых чисел — Решето Аткина

В процессе реализации одной программы я столкнулся с задачей поиска простых чисел до числа N порядка 10^9. На хабре уже неоднократно писали про различные способы и методы, но не упоминали про основной метод решения, что я сегодня и постараюсь исправить.



Алгоритм имеет асимптотическую сложность O(N/ln(ln(N)) и требует бит памяти. На входных значениях порядка 10^9 решето Аткина быстрее решета Эрастофена в 9.2 раза. Приведу график роста превосходства алгоритма Аткина на числах от 2 до 10^9:

График

В результате можно наблюдать следующую скорость выполнения:
10 000 000 0.15 сек.
100 000 000 2.16 сек.
1 000 000 000 48.76 сек.


Для удобства и краткости создадим абстрактный класс, от которого будем наследовать реализации решета Аткина. В дальнейшем будем только создавать различные варианты класса.
  public abstract class AAtkin
  {
    internal readonly int _limit;
    public bool[] IsPrimes;

    protected AAtkin(int limit)
    {
      _limit = limit;
      FindPrimes();
    }

    public abstract void FindPrimes();
  }



Создадим базовую реализацию алгоритма:
public override void FindPrimes()
    {
      IsPrimes = new bool[_limit + 1];
      double sqrt = Math.Sqrt(_limit);
        var limit = (ulong) _limit;
        for (ulong x = 1; x <= sqrt; x++)
          for (ulong y = 1; y <= sqrt; y++)
          {
            ulong x2 = x*x;
            ulong y2 = y*y;
            ulong n = 4*x2 + y2;
            if (n <= limit && (n%12 == 1 || n%12 == 5))
              IsPrimes[n] ^= true;

            n -= x2;
            if (n <= limit && n%12 == 7)
              IsPrimes[n] ^= true;

            n -= 2 * y2;
            if (x > y && n <= limit && n%12 == 11)
              IsPrimes[n] ^= true;
          }

        for (ulong n = 5; n <= sqrt; n += 2)
          if (IsPrimes[n])
          {
            ulong s = n*n;
            for (ulong k = s; k <= limit; k += s)
              IsPrimes[k] = false;
          }
      IsPrimes[2] = true;
      IsPrimes[3] = true;
    }



Время исполнения базовой (классической) версии алгоритма:
1 000 000 0.03 сек.
10 000 000 0.39 сек.
100 000 000 4.6 сек.
1 000 000 000 48.92 сек.


Вспомним, что в языке C# операции с типом ulong занимают примерно в 3 раза больше времени, чем с типом int. В программе при расчёте больших чисел значения всё-таки переваливают за int.MaxValue, так что найдём эмпирическую планку использования. Здесь она равна 858599509. Вставим условие и получим финальное ускорение:
IsPrimes = new bool[_limit + 1];
      double sqrt = Math.Sqrt(_limit);
      if (sqrt < 29301)
      {
        for (int x = 1; x <= sqrt; x++)
          for (int y = 1; y <= sqrt; y++)
          {
            int x2 = x * x;
            int y2 = y * y;
            int n = 4 * x2 + y2;
            if (n <= _limit && (n % 12 == 1 || n % 12 == 5))
              IsPrimes[n] ^= true;

            n -= x2;
            if (n <= _limit && n % 12 == 7)
              IsPrimes[n] ^= true;

            n -= 2 * y2;
            if (x > y && n <= _limit && n % 12 == 11)
              IsPrimes[n] ^= true;
          }

        for (int n = 5; n <= sqrt; n += 2)
          if (IsPrimes[n])
          {
            int s = n * n;
            for (int k = s; k <= _limit; k += s)
              IsPrimes[k] = false;
          }
      }
      else
      {
        var limit = (ulong) _limit;
        for (ulong x = 1; x <= sqrt; x++)
          for (ulong y = 1; y <= sqrt; y++)
          {
            ulong x2 = x*x;
            ulong y2 = y*y;
            ulong n = 4*x2 + y2;
            if (n <= limit && (n%12 == 1 || n%12 == 5))
              IsPrimes[n] ^= true;

            n -= x2;
            if (n <= limit && n%12 == 7)
              IsPrimes[n] ^= true;

            n -= 2 * y2;
            if (x > y && n <= limit && n%12 == 11)
              IsPrimes[n] ^= true;
          }

        for (ulong n = 5; n <= sqrt; n += 2)
          if (IsPrimes[n])
          {
            ulong s = n*n;
            for (ulong k = s; k <= limit; k += s)
              IsPrimes[k] = false;
          }
      }
      IsPrimes[2] = true;
      IsPrimes[3] = true;
    }



На этом всё. Спасибо за внимание — я постарался не быть очень сухим.
Поделиться публикацией

Комментарии 13

    0
    А где бывает прикладное применении этой задачи?
      +2
      Современная криптография использует именно этот алгоритм (с различными модификациями) для генерации последовательных простых чисел.
        0
        ясно, спасибо.
          0
          Ну это все понятно, куда интересней другой вопрос — а лично вы в каких целях использовали этот алгоритм?
            0
            Лично я — в учебных, если быть честным.
              +4
              В целях попадания на хабр, ясно же :-)
                0
                На GCJ один раз решето Аткина меня спасло.
                  0
                  Хм. На GCJ можно было бы как раз и не париться с генерацией простых чисел, а просто воспользоваться каким-нибудь сайтом с кучами этих чисел, слить оттуда многомегабайтный файлик, и воспользоваться им.

                  Какого года и какая именно задачка, если не секрет?
                    0
                    Использовать этот архив неспортивно:)
                    Не помню, хоть убейте. Помню, что написал его в template, где-то использовал, после контеста проверил — с решетом Эратосфена мог бы не успеть.
              0
              *применение
              +6
              Реализация есть — описания алгоритма нету.
              Минус Вам, уважаемый. :-)
              –1
              На вики алгоритм изложен короче и яснее, а приведенный там код на C дал для значения 8 000 000 с выводом всех чисел файл время 0m0.280s. Дальше правда не считает, наверно связано с тем, что использует тип int. Время расчета без вывода в файл 0m0.190s.

              Машина
              Architecture: x86_64
              CPU op-mode(s): 32-bit, 64-bit
              CPU(s): 2
              Vendor ID: GenuineIntel
              CPU family: 6
              CPU MHz: 2500.000
              L1d cache: 32K
              L1i cache: 32K
              L2 cache: 2048K

              Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

              Самое читаемое