Algorithm for finding the N first prime numbers - Atkin Sieve

Published on August 04, 2011

Algorithm for finding the N first prime numbers - Atkin Sieve

In the process of implementing one program, I was faced with the task of finding primes up to the number N of order 10 ^ 9. On a habr already repeatedly wrote about various ways and methods, but did not mention about the main method of the decision that I will try to correct today.



The algorithm has asymptotic complexity O (N / ln (ln (N))and requires a bit of memory. At the input values ​​of the order, the 10 ^ 9Atkin sieve is 9.2 times faster than the Erastofen sieve. Here is a graph of the growth of the superiority of the Atkin algorithm on numbers from 2 to 10 ^ 9:

Schedule

As a result, you can observe the following execution speed:
10,000,000 0.15 sec
100,000,000 2.16 sec
1,000,000,000 48.76 sec


For convenience and brevity, we will create an abstract class from which we will inherit the Atkin sieve implementations. In the future, we will only create various class variants.
  public abstract class AAtkin
  {
    internal readonly int _limit;
    public bool[] IsPrimes;

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

    public abstract void FindPrimes();
  }


Let's create a basic implementation of the algorithm:
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;
    }


The execution time of the basic (classic) version of the algorithm:
1,000,000 0.03 sec
10,000,000 0.39 sec
100,000,000 4.6 sec
1,000,000,000 48.92 seconds


Recall that in C #, operations with the ulong type take about 3 times longer than with the int type . In the program, when calculating large numbers, the values ​​still exceed int.MaxValue, so we will find an empirical usage bar. Here it is 858599509. Insert the condition and get the final acceleration:
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;
    }


That's all. Thank you for your attention - I tried not to be very dry.