I THINK ∴ I'M DANGEROUS

Prime Finder

I can't really remember why I decided I needed to write a program to find primes. But I did.

I used the GMP library so it will theoretically continue to run and find primes until it runs out of physical memory (rather than hitting some limitation of 64-bit integers). I am using an unsigned long long int for the index of the prime list. So that does limit me to 264 − 1 primes. If I actually hit this limit I may change to using a linked list.

Phases

  1. Highspeed, effiecient, Single threaded prime finder COMPLETE
  2. multithreaded primefinder
  3. distributed primefinder

Notes

Phase 1 was easy enough to complete. There may be some minor efficiency improvements to be made, but I'm basically done with it. Phase 2 has proven to be quite a challenge for me. I really like the idea of continuing to use the same basic algorithm as used in phase 1 (generate primes, use list of primes to quickly and deterministically identify prime numbers), however keeping the list sync'd and ordered across multiple threads has proven beyond my ability (and patience) to code.

After reading and digesting the wikipedia entry for Primality Tests, I am going to prototype my code in Tcl then port the algorithm into threaded C.

Code

Phase 1 code

primes.c
/*  Compile with gcc -o primes primes.c -O3 -lgmp
 */
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <gmp.h>
 
 
int main(void)
{
	unsigned long long p_index, p_max = 1;
 
	mpz_t p, p_sqrt;
	mpz_t *primes = malloc (sizeof(mpz_t) * p_max);
 
	mpz_init (p);
	mpz_init (p_sqrt);
	mpz_init (primes[0]);
 
	mpz_set_ui (p, 3);
	mpz_set_ui (primes[0], 3);
 
	for(;;)
	{
		next_p:
		mpz_sqrt (p_sqrt, p);
		mpz_add_ui (p_sqrt, p_sqrt, 1); /* perhaps unnecessary */
 
		for (p_index = 0; p_index < p_max; p_index++)
		{
			/* a number's factor will never be larger than its square root */
			if ( mpz_cmp (primes[p_index], p_sqrt) > 0 )
				break;
 
			if ( mpz_divisible_p (p, primes[p_index]) )
			{
				/* the number p is not prime */
 
				/* increment p by 2 */
				mpz_add_ui (p, p, 2);
 
				/* jump! */
				goto next_p;
			}	
		}
 
		/* if we got here that means p was not divisible by any 
		 * of the primes we have in our list which means p must
		 * be prime as well. Add it to the list */
 
		primes = realloc(primes, sizeof(mpz_t) * ++p_max);
		mpz_init (primes[p_max-1]);
		mpz_set (primes[p_max-1],p);
		mpz_add_ui (p, p, 2);
 
		/* show the number */
		mpz_out_str (stdout, 10, primes[p_max-1]);
		putchar (10);
	}
}