primes

kelso@gwusun.gwu.edu (John Kelso) (06/06/90)

There is a program on our Sun called primes which generates prime numbers.
I've seen it around on Unix systems for years.  Does anyone out there know
the algorithm that it uses, or better yet, where I can get the source for
it?  I'd read the fine manual, but we don't seem to have one for this
program.

Thanks,

      John Kelso, System Engineer, George Washington University
     SEAS Computing Facility, 725 23rd St NW, Washington DC 20052
 	     kelso@seas.gwu.edu  -or-  uunet!gwusun!kelso

falk@peregrine.eng.sun.com (Ed Falk) (06/07/90)

In article <8584@brazos.Rice.edu> kelso@gwusun.gwu.edu (John Kelso) writes:

|Does anyone out there know the algorithm that it uses, or better yet,
|where I can get the source for it?

It uses the "Sieve of Eratosthenes".  It's a very simple algorithm:

    1) make a list of N numbers to test:
	2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
    2) walk through the list, crossing out all multiples of two:
	2 3   5   7   9    11    13    15    17    19    21    23    25
    3) walk through the list, crossing out all multiples of three:
	2 3   5   7        11    13          17    19          23    25
    4) Do nothing for four, because it's already been crossed out.
    5) walk through the list, crossing out all multiples of five:
	2 3   5   7        11    13          17    19          23

    keep doing this for all n <= sqrt(N) that are still in the list when
    you get to them.

In going through my toolbox, I see that I have two different
implementations of this.  Apparently it was easier to just write it again
than remember what I called it the first time.

Note that you don't actually need to allocate enough room for a list of
numbers, just for a list of flags indicating that this number is crossed
out or not.

Another optimization is to simply ignore even numbers everywhere.

/* simple sieve of Eratosthenes for primes */

#include <stdio.h>
#include <math.h>

main(argc,argv)
	int	argc ;
	char	*argv[] ;
{
register unsigned char	*array ;
register int	i ;
	int	lim, n, j ;

	n = atoi(argv[1]) ;

	array = (unsigned char *) malloc(n+1) ;

	for(i=1;i<=n;++i)
	  array[i] = 1 ;

	lim = (int) sqrt((float) n) ;

	for(j=2; j<=lim; ++j)
	  if( array[j] != 0 )
	  {
	    for(i=2*j; i<=n; i += j)
	      array[i] = 0 ;
	  }

	for(i=2; i<=n; ++i)
	  if( array[i] != 0 )
	    printf("%d\n",i) ;
}
	-ed falk, sun microsystems -- sun!falk, falk@sun.com
	"What are politicians going to tell people when the
	Constitution is gone and we still have a drug problem?"
			-- William Simpson, A.C.L.U.