XGC - eXtended Eratosthenes Sieve

The Eratosthenes Sieve

We know the Eratosthenes Sieve as a method for selecting primes from the sequence of natural numbers. The algorithm of the sieve is very simple:

ES (for prime numbers): Initialize: you write down the natural sequence starting from 2

Sieve: you repeatedly scan the base sequence proceeding from less to greater indexes, select the first number p you find unmarked and mark a number each p

Terminate: at any moment you interrupt the scan, you can select all the unmarked numbers less than the square t of the least of them. All the selected numbers are all the prime numbers up to t

As all the primes except 2 are odd, we would like to find an Eratosthenes-like Sieve for applying only to odd numbers. Such a sieve not only exists, but it is still the Eratosthenes Sieve, provided that we change the "Initialize" step in this way:

ES for odd prime numbers: Initialize: you write down the 1 + 2 k sequence starting from 3

Sieve: you repeatedly scan the base sequence proceeding from less to greater indexes, select the first number p you find unmarked and mark a number each p

Terminate: at any moment you interrupt the scan, you can select all the unmarked numbers less than the square t of the least of them. All the selected numbers are all the prime numbers up to t

The ES practically didn't change because of the fact that each odd natural, that happens to be a multiple of a given odd prime p, is located (off from the previous multiple of p) at 2 p positions counting in the natural sequence, but only at p positions counting in the odd sub-sequence.
So, marking in the natural sequence a number each p corresponds to pre-marking all the numbers of the 2 k sub-sequence and marking in the 1 + 2 k sub-sequence a number each p.

Example

0: 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, ...
1: 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, ...

The line 0 lists all the even numbers and the line 1 all the odd ones. Selected primes are bold and marked multiples are underlined. The number 3 is the given odd prime, and the numbers 9, 15, 21, ... are the multiples of 3. As you can see,

9 is at 2 3 = 6 positions off from 3 if you count both lines
but it is at only 3 positions off if you count the line 1 only
15 too is at 6 positions off from 9 if you count both lines
but it is at only 3 positions off if you count the line 1 only
21 too is at 6 positions off from 15 if you count both lines
but it is at only 3 positions off if you count the line 1 only
...

We have to note that the "Sieve" step is highly optimized for the structure of the sequence given in the "Initialize" step, ie it is strictly bound to the concept of equally spaced terms in the base sequence.

There is a method for selecting exactly the same odd primes that doesn't use that feature. Instead, the alternate method use the concept of coprimality, ie the fact that two numbers have no common factors except 1. Of course, this concept is closely related to the definition of a prime number.
Testing two numbers for coprimality is not a difficult task and is far cheaper (in terms of computing resources) than testing a number for "simple" primality. In fact, testing for coprimality is as simple as calculating the Greatest Common Divisor (GCD) and, thanks to the Euclidean algorithm, the GCD is a simple task.
If the two numbers are a and b, then a and b are coprimes iif GCD[ a, b ] == 1

ES for odd prime numbers using coprimality: Initialize: you write down the 1 + 2 k sequence starting from 3

Sieve: you repeatedly scan the base sequence proceeding from less to greater indexes, select the first number p you find unmarked and mark each number n such that Not[ IsCoprime[ p, n ] ]

Terminate: at any moment you interrupt the scan, you can select all the unmarked numbers less than the square t of the least of them. All the selected numbers are all the prime numbers up to t

The eXtended Eratosthenes Sieve

The eXtended Eratosthenes Sieve is a method for selecting numbers from a given set of natural numbers. The starting set can be of any type, with or without a rule for determining any number of the set. By definition:

XES (general): Initialize: you write down the set

Sieve: you repeatedly scan the base set proceeding from less to greater indexes, select the first number p you find unmarked and mark each number n such that Not[ IsCoprime[ p, n ] ]

Terminate: all the selected numbers are all the coprime numbers of the set, scanning bottom-up.

Here is the formulation as a Mathematica function:

XES[ set ]

XES[ set_List ] :=
  Module[ {thisTerm, selectedTerms = {}, termsToSelect = set}, 
    While[ termsToSelect != {},
      thisTerm = First[ termsToSelect ];
      selectedTerms = {selectedTerms, thisTerm};
      termsToSelect = DeleteCases[ Rest[ termsToSelect ], eachTerm_ /; Not[ IsCoprime[ thisTerm, eachTerm ] ] ]
    ];
    Flatten[ selectedTerms ]
  ] /; IsNatural[ set ]

It gives the set resulting from the application of the eXtended Eratosthenes Sieve to the natural set.

Example

naturals = Range[ 2, 50 ]
{2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50}

primes = XES[ naturals ]
{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47}
oddNaturals = Range[ 3, 50, 2 ]
{3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 37, 39, 41, 43, 45, 47, 49}

oddPrimes = XES[ oddNaturals ]
{3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47}
congruence[ 2, 5 ] = Range[ 7, 100, 5 ]
{7, 12, 17, 22, 27, 32, 37, 42, 47, 52, 57, 62, 67, 72, 77, 82, 87, 92, 97}

XES[ congruence[ 2, 5 ] ]
{7, 12, 17, 37, 47, 67, 97}
noRule = Table[ Random[ Integer, {2, 100} ], {50} ]
{2, 18, 80, 33, 63, 72, 35, 5, 13, 92, 81, 3, 51, 5, 66, 24, 2, 52, 17, 92, 63, 2, 64, 74, 46, 20, 78, 14, 74, 95, 83, 37, 83, 86, 52, 50, 83, 89, 32, 38, 79, 77, 13, 66, 7, 86, 76, 71, 42, 96}

XES[ noRule ]
{2, 33, 35, 13, 17, 83, 37, 89, 79, 71}

If we know that the set is ordered and equally spaced then it is possible to optimize XES, just like ES, because of a nice property of that kind of sets.
Let's name m the space between adjacent terms. If the n-th term of the set is equal to f g, then each (n + f k)-th term is equal to f g + f k m = f (g + k m), ie the "Sieve" step has to mark each (n + f k)-th term of the set, as well as each (n + g k)-th one.

XES for equally spaced sets: Initialize: you write down the set

Sieve: you repeatedly scan the base set proceeding from less to greater indexes, select the first number p you find unmarked and, for each prime factor f of p, mark a number each f

Terminate: all the selected numbers are all the coprime numbers of the set, scanning bottom-up.

Here is the formulation as a Mathematica function:

XES[ set ]

pXES[ set_List ] :=
  Module[ {thisTerm, selectedTerms = {}, termsToSelect, primeFactors, thisFactor, len}, 
    len = Length[ set ];
    termsToSelect = Table[ True, {len} ];
    For[ t = 1, t <= len, t++,
      If[ termsToSelect[[ t ]],
        thisTerm = set[[ t ]];
        selectedTerms = {selectedTerms, thisTerm};
        primeFactors = First[ Transpose[ FactorInteger[ thisTerm ] ] ];
        While[ primeFactors != {},
          thisFactor = First[ primeFactors ];
          For[ i = t + thisFactor, i <= len, i += thisFactor,
            termsToSelect[[ i ]] = False
          ];
          primeFactors = Rest[ primeFactors ]
        ]
      ]    
    ];
    Flatten[ selectedTerms ]
  ]

XES[ set_List ] := pXES[ set ]  /; IsNatural[ set ] && OrderedQ[ set ] && IsEquallySpaced[ set ]

It gives the set resulting from the application of the eXtended Eratosthenes Sieve to the natural set.