/*	prime.c
	John Simpson <jms1@jms1.net> 2003-02-02

	quick-n-dirty prime number finder - uses 32-bit integers
	does not use any arrays
	does not use any floating point operations
	only checks each number for divisibility by previously-known primes

	this idea by Steve Litt <slitt@troubleshooters.com>.

	2003-02-04 - adding an optimization which pre-loads the first two
	primes (2 and 3) and then uses "add 2, add 4, add 2, add 4, ..."
	instead of just incrementing the number we're checking.

	2005-04-09 jms1 - changed copyright notice to indicate my intention   
	  that this code is licensed under GPL version 2 only, rather than 
	  GPL version 2 "or later". 

	****************************************

	Copyright (C) 2004-2005 John Simpson.

	This program is free software; you can redistribute it and/or modify
	it under the terms of the GNU General Public License, version 2, as
	published by the Free Software Foundation.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program; if not, write to the Free Software
	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
	or visit http://www.gnu.org/licenses/gpl.txt
*/

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>

#define NUM unsigned long
#define FORMAT "%lu"

NUM *plist ;
NUM *pend ;
int pcount ;

int main ( int argc , char *argv[] )
{
	int fcount , prime ;
	NUM n , r , s , t ;
	NUM *q ;
	unsigned short j ;

	struct timeval  ts , te ;
	struct timezone tz ;

	if ( argc < 2 )
	{
		puts ( "Count not specified" ) ;
		return 1 ;
	}

	fcount = atoi ( argv[1] ) ;

	if ( fcount < 1 )
	{
		puts ( "Count must be 1 or higher" ) ;
		return 1 ;
	}

	/* grab enough memory to hold all the primes we're gonna find */

	plist = (NUM *) malloc ( fcount * sizeof ( NUM ) ) ;
	if ( ! plist )
	{
		perror ( "malloc" ) ;
		return 1 ;
	}

	/* preload primes */

	pend = plist ;
	pcount = 0 ;

	*pend++ = 2UL ;
	pcount ++ ;

	if ( fcount > 1 )
	{
		*pend++ = 3UL ;
		pcount ++ ;
	}

	/* start the timer */

	gettimeofday ( &ts , &tz ) ;

	/* let's do it

		pcount	how many primes we've found
		plist	pointer to beginning of that list
		pend	pointer to end of that list

		n	the number we're testing for prime-ness
		r	integer "following" n
		s	the square of r (makes comparisons faster)
		t	already-known prime that we're testing against
		j	how far to jump to find the next number to test
			cycles between 2 and 4, which skips multiples
			of two and multiples of three.
	*/

	r = s = 1 ;
	j = 4 ;

	for ( n = 5 ; pcount < fcount ; n += ( j ^= 6 ) )
	{
		/*
			if n >= s, we know that sqrt(n) >= r
				need to increment r until s > n
		*/

		while ( n >= s )
		{
			r ++ ;
			s = r * r ;
		}

		/* check number against all previously-found primes */

		prime = 1 ;
		for ( q = plist ; prime && ( q < pend ) ; q ++ )
		{
			/* get the previously-known prime */

			t = *q ;

			/* if t > r, we know t > sqrt(n), therefore n is prime */

			if ( t > r )
			{
				break ;
			}

			/* if n is evenly divisible by t, n is not prime */

			if ( ! ( n % t ) )
			{
				prime = 0 ;
				break ;
			}
		}

		/* search is done. is n prime? */

		if ( prime )
		{
			/* add it to the end of the list */
			*pend++ = n ;

			/* increment the counter as well */
			pcount ++ ;

			/* maybe tell the user we found one */
			/* printf ( "(%d) " FORMAT "\n" , pcount , n ) ; */
		}
	}

	/* finish the timer */
	gettimeofday ( &te , &tz ) ;

	/* show the user what we found */
/*
	for ( q = plist ; q < pend ; q++ )
	{
		printf ( FORMAT "\n" , *q ) ;
	}
*/

	/* or maybe show them just the last prime */

	printf ( "Prime #%d is " FORMAT "\n" , pcount , plist[pcount-1] ) ;

	/* and show them how long it took */

	unsigned long long stim , etim , rtim ;
	stim = ts.tv_sec * 1000000 + ts.tv_usec ;
	etim = te.tv_sec * 1000000 + te.tv_usec ;
	rtim = etim - stim ;

	printf ( "Elapsed: %llu.%06llu sec\n" ,
		( rtim / 1000000 ) , ( rtim % 1000000 ) ) ;

	return 0 ;
}