Today is the birthday of Bernhard Riemann, a German mathematician who made fundamental contributions to the fields of geometry, analysis, and number theory. Riemann is definitely on my list of the greatest mathematicians of all time, and his conjecture about the distribution of prime numbers is one of the great unsolved problems in mathematics.
In honor of Riemann's birthday, this post deals with prime numbers.
One of the first mathematical algorithms that I learned as a kid is the “prime number sieve." (There is an animation of this algorithm in the Wikipedia article.) This algorithm is attributed to Eratosthenes, a mathematician from ancient Greece, and has been implemented in many programming languages. It is the only algorithm that I know of that has a poem written about it:
Sift the twos and sift the threes...,
The Sieve of Eratosthenes.
When the multiples sublime,
The numbers that remain are prime.
In celebration of Riemann’s birthday, I present the Sieve of Eratosthenes in the SAS/IML language:
/** The Sieve of Eratosthenes **/ proc iml; max = 200; list = 2:max; /** find prime numbers in [2, max] **/ primes = j(1, ncol(list), .); /** allocate space to store primes **/ numPrimes = 0; p = list[1]; /** 2 is the first prime **/ do while (p##2 <= max); /** search through sqrt(max) **/ idx = loc( mod(list, p)=0 );/** find all numbers divisible by p **/ list = remove(list, idx); /** remove them. They are not prime. **/ numPrimes = numPrimes + 1; primes[numPrimes] = p; /** include p in the list of primes **/ p = list[1]; /** next prime number to sift **/ end; /** include remaining numbers greater than sqrt(max) **/ k = numPrimes; n = k + ncol(list); primes[k+1:n] = list; /** put them at the end of the list **/ primes = primes[ ,1:n]; /** get rid of excess storage space **/ |
The following table shows the prime numbers less than two hundred that are found by the algorithm:
|
primes |
| 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, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199 |
Do YOU have a classic algorithm or technique that you have implemented in the SAS/IML language? (Poems are optional.)
4 Comments
proc iml;
**//generate the first N prime number//;
start prime(n);
A=2;
x=3;
k=1;
do until(k>=n);
do i=3 to x by 2 until(mod(x,i)=0);
if i=x then
do;
A=A//x; **//A collects all the n primes//;
k=k+1; **//k counts the number of primes;
end;
end;
x=x+1;
end;
return(A);
finish;
primes_n=prime(100);
print primes_n;
quit;
Pingback: Vectorized computations and the birthday matching problem - The DO Loop
Hi Rick, The following is amazingly inefficient but it can be done in one-line of IML! If the task is to find only a small number of primes, then why not multiply everything together (twice) and compare to the list of postive integers. It might even run faster than the sieve, but only if you have a max n of about 50 or less!
proc iml;
start primes_upto(n);
return(setdif( 2:n, t(2:int(n/2))*(2:int(n/2)) ));
finish;
print (primes_upto(200));
quit;
I created my own prime number sieve.
This is a different prime number sieve. It omits all the even numbers and eliminates all the odd numbers which are non-prime. Here's how the sieve works. We begin with the first odd number after the number 1. That is 3. We perform 1 addition operation at the number 3. The addition operation is as follows. We take the number immediately preceding 3 which 2 and add it to 3. Then we take the number immediately succeeding 3 which is 4 and add it to our total like so, 2+3+4=9. The purpose of this operation is to tell us the fist odd number in sequence we will arrive at after 3 which will not be a prime number. In other words, the odd numbers in between 3 and 9(those are 5 and 7) are prime numbers while the number 9 is not. Not only that if we square root our result,the number 9, we also find out that we arrive at the first prime number after 2 which is the sqaure root of 9 which is 3. In other words, our operation here at prime number 3 tells which numbers will not be prime and if we square root our result it confirms our first prime number after 2.
Now that we have performed these operations at the first odd number 3, we move on to the next odd number in sequence after 3 which is 5(we already know it's a prime number but we have to perform the operations here at 5 anyways). Only this time we are going to do not one but two mathematical additions which follow this simple pattern. Just like with 3 our first operation takes the number immediately precededing 5 which is 4 and the number immediately succeeding 5 which is 6 and adds them to 5 like so 4+5+6=15. 15 is the next odd number we will encounter in our succession of odd numbers that will not be a prime number. In other words, all odd numbers from 3 to 15 with the exception of 9 and 15 are prime numbers. Now that we have performed one operation at 5 we perform a second one that follows the same pattern. Our first operation was to take the number immmmediately preceding 5 and immediately succeeding 5 and add the total. This next operation is the same only we go one step further. We now take the TWO numbers immediately preceding 5 which is 3 and 4 and the TWO numbers immediately succeeding 5 which 6 and 7 and we add all of them like so, 3+4+5+6+7=25. 25 will also not result in a prime number and also notice this, if we square the result we end up with 5 again.
Now that we have performed these operations at 3 and 5 and have eliminated the following odds numbers as not prime(9,15,25), we move onto the next odd number in sequence which is 7. Care to guess how many mathematical operations we perform here? Well, we continue with the pattern and perform 3 of them. Our first operation is 6+7+8=21. 21 will will not be a prime number. The next operation is 5+6+7+8+9=35. 35 will not be a prime number. The third and final operation at 7 will be 4+5+6+7+8+9+10=49(where we take the 3 preceding numbers and the three succeeding numbers). 49 will not be a prime number and if we sqaure the result we find that the square root of 49 is 7.
Now let's continue this sieve to show how it continues to work. The next odd number after 7 is 9. We've already previously learned that 9 is not a prime number but we can't skip the operations at it for the purposes of the sieve. We need to perform 4 operations at 9 and every succeeding odd number thereafter will increment the number of operations by 1 each time. So our first operation here will be, 8+9+10=27. 27 will not be a prime number and is eliminated from our odd number pool of possible primes as with our previous operations. The next operation is 7+8+9+10+11=45. 45 is not a prime number. The next operation is 6+7+8+9+10+11+12=63. 63 is not a prime number. And the final operation here is 5+6+7+8+9+10+11+12+13=81. 81 is not a prime. Since we already know that 9 is not a prime number from previous operations there is no need to squareroot 81 here. We move onto the next odd number after 9 which is 11 and perform 5 operations there and on and on through the odd number progression.
That is how the sive works and I believe we may be able to quantify this sieve with a simple equation but how to do that I am unsure of.