Sundarama Sieve

Sundarama sieve in the network is represented by a large number of sources of brief background information. Nevertheless, I decided to state what I myself would like to read at the beginning of the study of number-theoretic algorithms.



Sundarama sieve is one of the three most famous methods for generating prime numbers. Now it is customary to treat it as some exotic because of poor computational complexity: O (N (logN)). However, asymptotics are asymptotics, and in practice, in the 32-bit sieving range, Atkin, for example, overtakes Sundaram only with careful optimization.



Implementations of the Atkin sieve, circulating on the Internet, do not surpass Sundaram sieve neither in terms of time, nor in memory efficiency. So the Sundaram method can be used as an auxiliary tool in experiments with more advanced algorithms.



Sundarama sieve finds all primes in a given range of natural numbers 3 ≀ n ≀ N "sifting out" the components. Without loss of generality, the value of N can be considered odd. If N is even, then it is guaranteed to be composite, and it can be excluded from the sieving range by decreasing the value of the upper boundary by one.



The algorithm is based on the following fact. Any odd composite number n can be represented as the product of two natural odd numbers greater than one:



(1) n = (2 * i + 1) * (2 * j + 1),



where i and j are natural numbers (zero is not a natural number). It is impossible to imagine a prime number n in this form, since otherwise it would mean that n has additional divisors besides himself and one.



We write the odd n in the form 2 * m + 1, substitute it in (1), and get the expression for m:



(2) m = 2 * i * j + i + j.



This transformation directly leads to the idea of ​​the sieve of Sundaram.



In order to filter out composite numbers from a given interval, you should use an array in which an element with index m is a representative of the number 2 * m + 1. Having organized the enumeration of the values ​​of the variables i and j, we will calculate the index m by the formula (2), and in the corresponding array elements set the mark "composite number". After the enumeration is completed, all the composite numbers in the range will be marked, and primes can be selected by the absence of a mark.



It remains to clarify the technical details.



Based on the previous reasoning, to represent the upper (odd) boundary of the sieving range N, the index m assumes its maximum value mmax = (N - 1) / 2. This determines the required size of the array.



We will enumerate the values ​​of i and j in two cycles: an outer loop for enumerating the values ​​of i, and a nested inner loop for the values ​​of j.



The initial value of the external loop counter is i = 1. Given the complete symmetry of the occurrence of the variables i and j in expression (2), to avoid repeated duplicate calculations, the internal cycle should begin with the value j = i.



Find the maximum value for the external loop counter imax β‰₯ i. If the boundary of the range N is an odd composite number, then with the value i = imax, the inner loop must be executed at least once with the value of its parameter j = imax to delete N, and expression (2) will take its maximum value:



mmax = 2 * imax * imax + imax + imax,

imax ^ 2 + imax - mmax / 2 = 0.



Solving this quadratic equation, we get: imax = (√ (2 * mmax + 1) -1) / 2 = (√N-1) / 2.

The restriction for the counter of the inner cycle jmax β‰₯ j is found from the inequality

mmax β‰₯ 2 * i * j + i + j , whence we get: jmax = (mmax - i) / (2 * i + 1).



The values ​​of the upper limits should be rounded to the nearest whole value, discarding the fractional part.



The following is a listing of a very simple C # Sundaram class that implements the described method.



using System; using System.Collections; namespace CSh_Sundaram { public class Sundaram { public double dt; //   () private long t1; //   private long t2; //   private uint limit; //     private int arrlength; //   private BitArray Prim; //    private int counter; public Sundaram(uint _limit) { limit = _limit; if (limit % 2 == 0) limit -= 1; arrlength = (int)(limit / 2); Prim = new BitArray(arrlength); t1 = DateTime.Now.Ticks; Sieve(); //  t2 = DateTime.Now.Ticks; dt = (double)(t2 - t1) / 10000000D; counter = -1; } public uint NextPrime { get { while (counter < arrlength - 1) { counter++; if (!Prim[counter]) return (uint)(2 * counter + 3); } return 0; } } private void Sieve() { int imax = ((int)Math.Sqrt(limit) - 1) / 2; for (int i = 1; i <= imax; i++) { int jmax = (arrlength - i) / (2 * i + 1); for (int j = i; j <= jmax; j++) { Prim[2 * i * j + i + j - 1] = true; } } } } }
      
      





As a parameter when creating an object of type Sundaram, the upper limit of the sieving range is indicated. For sifting, the class uses an array of type BitArray. This increases the runtime, but allows you to cover the entire 32-bit range of the uint type. Sifting is performed when accessing the Sieve () method from the class constructor. After its completion, the dt field will contain the Sieve execution time in seconds.



When accessing the NextPrime property, sequentially, starting from the value 3, in ascending order, primes will be returned. After all the simple ones from the range are exhausted, the value 0 will be returned.



All Articles