APL Wiki
The Sieve of Eratosthenes
Contents
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 0We 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 0Finally 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