Using Sieve of Eratosthenes to find prime numbers

Sep 25, 2017 18:56 · 499 words · 3 minutes read R Python Algorithm Sieve of Eratosthenes

Mathedemo

We know a prime number is a natural number greater than 1 has no positive divisors other than 1 and itself. If Someone asks you, “What are all prime numbers less than 1 million?”, how can you answer the question?

There are many algorithms to find all prime numbers up to any given limit. I’m going to cover the Sieve of Eratosthenes, a simple but powerful algorithm, in this article.

Sieve of Eratosthenes

To find all the prime numbers less than or equal to a given integer n by Eratosthenes’ method:

1.Create a list of consecutive integers from \(2\) through \(n\) \((2, 3, 4, ..., n)\).

2.Initially, let \(p = 2\), the smallest prime number.

3.Enumerate the multiples of \(p\) by counting to \(n\) from \(p ^ 2\) in increments of \(p\), and mark them in the list (these will be \(p ^ 2, p ^ 2 + p, p ^ 2 + 2p, ...\); the \(p\) itself should not be marked).

4.Find the first number greater than \(p\) in the list that is not marked. Let \(p\) now equal this new number (which is the next prime), and repeat from step 3. If \(p ^ 2\) is greater than \(n\), stop.

5.When the algorithm terminates, the numbers remaining not marked in the list are all the primes number not greater than \(n\).

Pseudocode

About algorithm can be expressed in pseudocode as below:

Input: an integer n > 1

Let A be an array of Boolean value, indexed by integers 2 to n, initially all set to TRUE.

for i = 2, 3, 4, ..., n ^ (1 / 2):
  if A[i] is TRUE:
    for j = i ^ 2, i ^ 2 + i, i ^ 2 + 2i, ..., n:
      A[j] = FALSE

Output: all i such that A[i] is TRUE.

Let’s see below gif illustrating the process:

Sieve of Eratosthenes: algorithm steps for primes below 121

Implementation in R

all_prime <- function(n) {
  prime <- rep(TRUE, n)  # Intitialize a vector all set to TRUE
  prime[1] <- FALSE  # Remove 1 from output
  
  for (i in 2:sqrt(n)) {
    if (prime[i]) prime[seq(i ^ 2, n, i)] <- FALSE  # Mark FALSE if composite index.
  }
  
  which(prime)  # Return indices of TRUE
}

# Example, find all prime numbers less or equal to 61
all_prime(61)

##  [1]  2  3  5  7 11 13 17 19 23 29 31 37 41 43 47 53 59 61

Implementation in Python

import math
import numpy as np

def all_prime(n):
    prime = np.array([True] * (n + 1))  # Intitialize a vector all set to TRUE
    prime[0] = prime[1] = False  # Remove 0 and 1 from output

    for i in range(2, int(math.ceil(math.sqrt(n)))):
        if prime[i]:
            prime[range(i ** 2, n, i)] = False  # Mark FALSE if composite index.

    print(np.where(prime == 1))  # Return indices of TRUE

all_prime(61)

## (array([ 2,  3,  5,  7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59,
##        61]),)