Finding primes from 1 upwards

[mdettweiler]

Originally Posted by henryzz

has anyone seached exhautively for primes particularly high
i have seached the web but cannot find any projects

if not what would be the best method for starting this search

The PrimeGrid BOINC project is already doing this with their "Primegen" subproject.

[ShiningArcanine]

Originally Posted by petrw1

I wrote a program a few years ago based on the Seive Of Eratosthenes. Due to memory limitations of the program I used (BASIC ) I had to stop at about the number 6 Billion. It only took about a day to run on a P3 400Mhz.

There are some links to lists here: http://primes.utm.edu/lists/small/

Recently, I computed all of the primes from 3 to 2^31-1 on my computer using a program I wrote that utilizes the sieve of erathothenes. My computer has a Core 2 Duo E6300 and 4GB of RAM (only 3.25 GB are accessible due to limitations of the operating system). According to my calculations, neglecting the memory requirements of the operating system, it should be possible for a modified version of the program I wrote to compute all of the primes between 1 and 223,270,238,621, with the exception of 2, which is not truly computed as I hard coded it in my program so it could neglect all of the even integers.

I posted my program on channel9.msdn.com in a thread I started to find out why it did not compile in Visual Studio (as I originally developed it using GCC on Unix). I subsequently posted later, more enhanced versions of it for additional criticism, as I wrote it to become more proficient in C and I wanted tips from more experienced programmers on how it could have been written to execute more quickly. I am at the university right now, so I do not have access to my copy of the program at home, but here is the most recent version I posted online:

/*  prime.c */
/*  Sept 29, 2007. */
/*  <my name was here> */
/*  Program to find prime numbers */
#include <stdlib.h>
 
int main (void)
{
 
    unsigned int len, b_num, num, bitVal, *arr = NULL, *arrEnd = NULL;
 
    printf ("\nFind prime numbers smaller than \n");
    scanf("%i", &len);
 
    if (len < 1)
    {
 
        return 0;
 
    }
 
    b_num = ( len % (16 * sizeof(unsigned int)) == 0 ) ?
        (len / (16 * sizeof(unsigned int))) :
        (len / (16 * sizeof(unsigned int)) + 1);
 
    arr = (unsigned int*)calloc (b_num, sizeof(unsigned int));
 
    arrEnd = arr + b_num;
 
    bin_arr_prime_fill(arr, b_num);
 
    if (len >= 2)
    {
 
        printf("Number 2 is a prime number\n");
 
    }
 
    /* Loop through Prime Array and Output Primes */
    for ( num = 1; arr < arrEnd ; ++arr )
    {
 
        for ( bitVal = 0x80000000; bitVal > 0 && num <= len; bitVal >>= 1, num += 2 )
        {
 
            if ((*arr & bitVal) == bitVal)
            {
 
                printf("Number %i is a prime number\n", num);
 
            }
 
        }
 
    }
 
    free(arr - b_num);
 
    main();
 
    return 0;
 
}
 
int bin_arr_prime_fill (const unsigned int * const arr, const unsigned int len)
{
 
    unsigned int b, num, *a = arr, *aEnd = arr + len;
 
    /* b has a binary 1 where the number is in the array */
    /* num is the number being checked for primiality */
    /* a is used to seek to the appropriate position in the array * /
    /* aEnd is used to determine when to stop */
    for ( num = 3, b = 0x40000000; a < aEnd; ++a, b = 0x80000000 )
    {
        for (  ; b > 0 ; b >>= 1, num += 2 )
        {
 
            /* arr is the prime array */
            /* a points to the block currently being built */
            /* num is the number being tested */
            if (isPrime(arr, num) == 1)
            {
                *a |= b;
            }
 
        }
 
    }
 
    return 0;
 
}
 
int isPrime(const unsigned int * arr, const unsigned int num)
{
 
    unsigned int bitVal = 0x40000000, div = 3;
 
    /* arr points to an array; it is used to seek to the correct block */
    /* num is the number being checked for primiality */
    /* bitVal is the bit position */
    /* div is the divisor, starting with 3 */
    /* Set bitVal = 0x2000 so it corresponds to m = 3 for the */
    /* first inner loop */
    /* Set bitVal = 0x8000 so it corresponds to the first number */
    /* in the block in every suceeding loop */
    /* The outer loop loops through each block of the array */
    for (  ; div <= (num / div); ++arr, bitVal = 0x80000000)
    { 
 
        /* The inner loop checks each number in the block for\ */
        /* primality, skipping even numbers */
        for (  ; div <= (num / div) ; bitVal >>= 1, div += 2 )
        {
 
            /* If the divisor is prime and it divides */
            /* n evenly, return composite */  
            if ((*arr & bitVal) == bitVal && (num % div) == 0)
            {
 
                return 0;
 
            }
 
        }
 
    }
 
    /* Return Prime */
    return 1;
 
}

Looking through it, it should be almost identical to the copy I have at home, which has a newer version main() function that uses fewer variables.

If my memory serves me correctly, it took 12 to 16 hours to do such a computation on my computer, although I did not take precise time measurements as I only wanted a ball park estimate of how long it took to compute. The only potential single threaded method I can think of doing to compute this more quickly would be to take advantage of differences in the computational speed of a processor's ALUs and FPUs by using floats instead of unsigned integers, but I have not tried it as I am working on another project. I am curious, how did you compute all of the primes between 1 and 6 billion in 24 hours on a 400MHz Pentium 3?

[axn]

Originally Posted by ShiningArcanine

If my memory serves me correctly, it took 12 to 16 hours to do such a computation on my computer

Yeesh! What you have coded is not Sieve of Erat. It is just brute force trial division with a bitmap to store the primes found so far.

Please take a look at here for an explanation of SoE.

A properly implemented SoE can do the range (2^31) in a few seconds (not hours).

[ShiningArcanine]

Originally Posted by axn1

Yeesh! What you have coded is not Sieve of Erat. It is just brute force trial division with a bitmap to store the primes found so far.

Please take a look at here for an explanation of SoE.

A properly implemented SoE can do the range (2^31) in a few seconds (not hours).

I have read that before today and I did not think that it made a difference whether or not I started with an bitmap of 0s and marked off the 1s or an bitmap of 1s and marked off the zeros. I guess because of the need to do integer division when marking off the 1s, it does. I will rewrite my program to initialize a bitmap with 1s and then mark off the 0s.

[ShiningArcanine]

I rewrote my program. It required far less code to start with a bitmap of ones and then mark off the zeroes. It also can sieve from 1 to 1073741824 in next to no time at all.

Here is the code:

/*  prime.c */
/*  Oct 17, 2007. */
/*  <my name was here> */
/*  Program to find prime numbers */
#include <stdlib.h>
	
int main (void)
{

	unsigned int len, b_num, num, bitVal, *arr = NULL;

	printf ("\nFind prime numbers smaller than \n");
	scanf("%i", &len);
			
	if (len < 1)
	{
				
		return 0;
			 
	}
		 
	b_num = ( len % (16 * sizeof(unsigned int)) == 0 ) ?
		(len / (16 * sizeof(unsigned int))) :
		(len / (16 * sizeof(unsigned int)) + 1);
	
	arr = (unsigned int*)calloc (b_num, sizeof(unsigned int));

	arr[0] = 0x7FFFFFFF;

	for ( num = 1 ; num < b_num ; ++num )
	{
	
		arr[num] = 0xFFFFFFFF;

	}
	
	bin_arr_prime_fill(arr, b_num);
 
	if (len >= 2)
	{

		printf("Number 2 is a prime number\n");

	}

	/* Loop through Prime Array and Output Primes */
	for ( num = 1; num <= len ; ++arr )
	{

		for ( bitVal = 0x80000000; bitVal > 0 && num <= len; bitVal >>= 1,
num += 2 )
		{
	
			if ((*arr & bitVal) == bitVal)
			{   
	 
				printf("Number %i is a prime number\n", num);
		
			}
	 
		}
	
	}
		
	free(arr - b_num);
	
	main();
	
	return 0;
	
}
	
int bin_arr_prime_fill (unsigned int * const arr, const unsigned int b_num)
{

	unsigned int num, i, len = b_num * 64;

	/* num is the number being checked for primiality */

	for ( num = 3; num * num <= len ; num += 2 )
	{

		if ((arr[num >> 6] & (0x80000000 >> ((num >> 1) % 32))) == 0) 
		{

			continue;

		}
		
		for ( i = num * 3; (i >> 6) < b_num ; i += num << 1 )
		{

			if (arr[i >> 6] & (0x80000000 >> ((i >> 1) % 32)))
			{
			
				arr[i >> 6] ^= (0x80000000 >> ((i >> 1) % 32));

			}

		}

	}
	
	return 0;
	
}

I had heard in the past that division operations were slow on computer processors, but the orders of magnitude difference I am seeing in performance here is ridiculous.

[Xyzzy]

http://www.troubleshooters.com/codec...imenumbers.htm

[bsquared]

Getting rid of the function call in the innermost loop certainly helped too. That was a fast rework of your code, good job.

[ShiningArcanine]

Originally Posted by bsquared

Getting rid of the function call in the innermost loop certainly helped too. That was a fast rework of your code, good job.

Thanks. I wrote this to get better at writing C code so reading that really means something to me.

[bsquared]

It looks like you got rid of all the division operators in the inner loops. You can also get rid of the % operators when %'ing by a power of 2. since division by a power of 2 (in this case 32) just shifts out the lowest bits, using & accomplishes the same thing. % 32 becomes & 31.

[bsquared]

I went ahead and tried out the above suggestion on your code because I'm looking for excuses to not do my homework at the moment

It had negligible impact (Compiling with MSVC 6).

To further increase speed (assuming you want to) you'll have to do something about the cache performance of the sieve. Read about cache blocking, or a segmented sieve.

This kind of optimization can be addictive. I've spent a fair amount of time tuning mine up, and while it's fast, it's still not the fastest I've seen (but close). There's always something else to try...

[ShiningArcanine]

I found a slight bug in my code, which lowered upper bound for the sieve fairly substantially. I had only tested very large upper bounds, which caused me to miss the bug. With the fix, for 1 to 1073741824, my program is not quite as fast as I thought it was earlier, but it is still substantially faster than what it was, taking 10 to 11 seconds to find all of the primes between 1 and 1073741824 on my computer versus the approximate 12 to 16 hours it took before the rewrite.

Thanks for the tip regarding using binary AND instead of the arithmetic modulus. I have made that change as well as a few other enhancements to my code. Among them I added a new function to make the code more readable. I have the /Ob2 flag set in Visual Studio instructing it to inline functions whenever it can, so I assume that it is inlined sometime during the compilation process. Here is the code for the new version:

/*  prime.c */
/*  Oct 17, 2007. */
/*  <my name was here> */
/*  Program to find prime numbers */
#include <stdlib.h>

int main (void)
{

	unsigned int len, num, bitVal, *arr = NULL;

	printf ("\nFind prime numbers smaller than \n");
	scanf("%i", &len);
			
	if (len < 1)
	{
				
		return 0;
			 
	}
	
	arr = (unsigned int*)calloc (blocks(len), sizeof(unsigned int));
	
	bin_arr_prime_fill(arr, len);
 
	if (len >= 2)
	{

		printf("Number 2 is a prime number\n");

	}

	/* Loop through Prime Array and Output Primes */
	for ( num = 1; num <= len ; ++arr )
	{

		for ( bitVal = 0x80000000; bitVal > 0 && num <= len; bitVal >>= 1,
num += 2 )
		{
	
			if ((*arr & bitVal) == bitVal)
			{   
	 
				printf("Number %i is a prime number\n", num);
		
			}
	 
		}
	
	}
		
	free(arr - blocks(len));
	
	main();
	
	return 0;
	
}

unsigned int blocks (const unsigned int len)
{

	return ( len & ((sizeof(unsigned int) << 4) - 1) == 0 ) ?
		(len / (sizeof(unsigned int) << 4)) :
		(len / (sizeof(unsigned int) << 4) + 1);

}

int bin_arr_prime_fill (unsigned int * const arr, const unsigned int len)
{

	unsigned int num, i, b_num = blocks(len);

	arr[0] = 0x7FFFFFFF;

	for ( num = 1 ; num < b_num ; ++num )
	{
	
		arr[num] = 0xFFFFFFFF;

	}

	/* num is the number being checked for primiality */

	for ( num = 3; num * num <= len ; num += 2 )
	{

		if ((arr[num >> 6] & (0x80000000 >> ((num >> 1) & 0x1F))) == 0) 
		{

			continue;

		}
		
		for ( i = num * 3; (i >> 6) < b_num ; i += num << 1 )
		{

			if (arr[i >> 6] & (0x80000000 >> ((i >> 1) & 0x1F)))
			{
			
				arr[i >> 6] ^= (0x80000000 >> ((i >> 1) & 0x1F));

			}

		}

	}
	
	return 0;
	
}