The Sieve of Eratosthenes

The Sieve in APL

The Sieve of Eratosthenes is a simple, well-known algorithm for finding Prime Numbers. It was created by the Ancient Greek mathemetician Eratosthenes.

The algorithm can be written in a single line of APL code which finds all the prime numbers up to and including X:

      (2=+⌿0=(⍳X)∘.|⍳X)/⍳X

For example, to list all the prime numbers up to 100:

      (2=+⌿0=(⍳X)∘.|⍳X)/⍳X←100
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97

The Sieve in the C programming language

It's interesting to compare the brevity of this single line of code with a typical implementation written in the C programming language. In C, the code might look as follows:

/* Sieve of Eratosthenes in C */
#include <stdio.h>
#include <stdlib.h>

int main (int argc, char **argv)
{
    unsigned long n, x, y, *primes;
 
    /* Get the upper limit value, n */   
    if (argc != 2) {
        fprintf (stderr, "Usage is e.g:\n    %s 10\n", argv[0]);
        return -1;
    }
    
    n = strtoul (argv[1], NULL, 0);
    if (n == 0) {
        fprintf (stderr, "Argument must be greater than 0\n");
        return -1;
    }
    
    /* Run the sieve algorithm */
    primes = (unsigned long *) calloc (n+1, sizeof (unsigned long));
    if (primes == NULL) {
        fprintf (stderr, "Out of memory\n");
        return -1;
    }
 
    for (x = 2; x <= n; x++) {
        for (y = 2; y <= x; y++) {
            if (x * y > n)
                break;
            primes [x * y] = 1;
        }
    }
    
    /* Print the results */
    for (x = 2; x <= n; x++) {
        if (primes [x] == 0)
            printf ("%d ", x);
    }
    printf ("\n");
    return 0;
}

And of course in C we would need to compile and link the program before we could run it. In APL, you just type in the expression.

How the APL version works

The APL version starts from the observation that if a number M is a multiple of another number N, then the remainder of M ÷ N must be zero. The APL function Residue | gives remainder values, e.g. 3|10 gives 1 so 10 is not evenly divisible by 3.

We can find the remainders when a whole vector of numbers is divided by another number. For example:

      3|1 2 3 4 5 6 7 8 9 10
1 2 0 1 2 0 1 2 0 1

And we can use APL's Outer Product ∘. together with Residue to find the remainder when the vector of numbers is divided by each number in another vector:

      3 4 5∘.|1 2 3 4 5 6 7 8 9 10      
1 2 0 1 2 0 1 2 0 1
1 2 3 0 1 2 3 0 1 2
1 2 3 4 0 1 2 3 4 0

So we can begin to understand the APL solution. Let's look for prime numbers up to 10:

      X←10
      ⍳X
1 2 3 4 5 6 7 8 9 10
      (⍳X)∘.|⍳X      
0 0 0 0 0 0 0 0 0 0
1 0 1 0 1 0 1 0 1 0
1 2 0 1 2 0 1 2 0 1
1 2 3 0 1 2 3 0 1 2
1 2 3 4 0 1 2 3 4 0
1 2 3 4 5 0 1 2 3 4
1 2 3 4 5 6 0 1 2 3
1 2 3 4 5 6 7 0 1 2
1 2 3 4 5 6 7 8 0 1
1 2 3 4 5 6 7 8 9 0

We now want to find out which divisions resulted in a zero remainder, which we can do by using the test 0 = . A '1' in the result for column M and row N indicates that M was divisible by N

      0=(⍳X)∘.|⍳X      
1 1 1 1 1 1 1 1 1 1
0 1 0 1 0 1 0 1 0 1
0 0 1 0 0 1 0 0 1 0
0 0 0 1 0 0 0 1 0 0
0 0 0 0 1 0 0 0 0 1
0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 1

A number M is prime if it's only divisible by itself and one. Such numbers will only have two 1s in column M. So the next thing to do is to sum the columns. Where the result is 2, the number is prime.

We can sum the numbers in each column by using the + function together with the First Axis Reduction operator, . This operator applies the + function along the first axis of the table (the rows) - i.e. for each column it adds all the row values together.

      +⌿0=(⍳X)∘.|⍳X
1 2 2 3 2 4 2 4 3 4
      2=+⌿0=(⍳X)∘.|⍳X
0 1 1 0 1 0 1 0 0 0

Finally we can take the vector of 0s and 1s and use APL's Compression function / to find the column positions corresponding to prime numbers:

      (2=+⌿0=(⍳X)∘.|⍳X)/⍳X
2 3 5 7

   

Author: SimonMarsden


CategoryShowcases

SieveOfEratosthenes (last edited 2010-01-08 19:52:56 by SimonMarsden)