Finding primes from 1 upwards

[axn]

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

can be replaced by

arr[i >> 6] &= ~(0x80000000 >> ((i >> 1) & 0x1F))

[axn]

Also, while clearing of multiples of num, instead of starting from num * 3, you can start from num * num (Ask yourself why this would work!).

Also, the whole inner loop can be optimised by directly working with j = i/2 rather than with i. Some thing like:

for ( i = num * num, j=i>>1; (j >> 5) < b_num; j += num)
{
    arr[j >> 5] &= ~(0x80000000 >> (j & 0x1F))
}

[axn]

Originally Posted by ShiningArcanine

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.

It is not just that the division operator is slow -- you can consider a div instruction to be 10x slower than a memory write, which would not explain the gains.

It is that, earlier, you were doing trial division -- i.e you were trying to divide every number with primes < sqrt(number), whereas SoE knows exactly which numbers are multiples of a given prime and checks off those only.

So, for example, your first version will try to check every single number to see if it is divisible by 3, whereas SoE will only check off the actual multiples of 3.

In fact, it'd be a very interesting exercise to count the inner loops of both versions -- see how many time SoE executes the "clear" operation, vs how many times the first version executed the "(num % div) == 0" operation.

Be warned, you'll probably need to use 64 bit ints to fit the count

[ShiningArcanine]

Originally Posted by axn1

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

can be replaced by

arr[i >> 6] &= ~(0x80000000 >> ((i >> 1) & 0x1F))

That does not work, as it sets all of the bits as equal to zero, except the bit being considered. However, your next suggestion will enable me, to use:

arr[i >> 6] ^= ~(0x80000000 >> ((i >> 1) & 0x1F))

Originally Posted by axn1

Also, while clearing of multiples of num, instead of starting from num * 3, you can start from num * num (Ask yourself why this would work!).

Also, the whole inner loop can be optimised by directly working with j = i/2 rather than with i. Some thing like:

for ( i = num * num, j=i>>1; (j >> 5) < b_num; j += num)
{
    arr[j >> 5] &= ~(0x80000000 >> (j & 0x1F))
}

My use of the following was actually a workaround because of this:

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

I think that j needs to be used slightly differently because this program does not have memory structures to test even integers.

So the inner loop becomes:

for ( i = num * num, j = num << 1; (i >> 6) < b_num ; i += j )
{

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

}

Originally Posted by axn1

It is not just that the division operator is slow -- you can consider a div instruction to be 10x slower than a memory write, which would not explain the gains.

It is that, earlier, you were doing trial division -- i.e you were trying to divide every number with primes < sqrt(number), whereas SoE knows exactly which numbers are multiples of a given prime and checks off those only.

So, for example, your first version will try to check every single number to see if it is divisible by 3, whereas SoE will only check off the actual multiples of 3.

In fact, it'd be a very interesting exercise to count the inner loops of both versions -- see how many time SoE executes the "clear" operation, vs how many times the first version executed the "(num % div) == 0" operation.

Be warned, you'll probably need to use 64 bit ints to fit the count

It would be fairly humorous if I used a 32bit int and due to the overflow, the numbers came back approximately the same.

By the way, here is the latest version of my code:

/*  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) == 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, j, 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 )
	{
		
		for ( i = num * num, j = num << 1; (i >> 6) < b_num ; i += j )
		{

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

		}

	}
	
	return 0;
	
}

I had to roll back part of the blocks() function because of a bug that was introduced when I removed the modulus.

[ShiningArcanine]

It seems that the inner loop's test condition can be changed from (i >> 6) < b_num to i < len, which eliminates the need for a bitshift every iteration of that loop.

[ShiningArcanine]

I had to bring back:

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

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

}

I had been wrong in thinking that the other optimization eliminated the need for the conditional, as 63 gets set to true for the range 1 to 64 because 63 is a multiple of 3 and 7.

Edit: With this, the inner loop becomes:

for ( i = num * num, j = num << 1 ; i < len ; i += j )
{
			
	if (arr[i >> 6] & (0x80000000 >> ((i >> 1) & 0x1F)))
	{

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

	}

}

Edit: I really should find time to comment this code. It seems I removed an optimization when I was trying to further optimize my code. Here is the new code:

/*  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) == 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, j, k, 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 * num, j = num << 1; i < len ; i += j )
		{
			
			if (arr[i >> 6] & (0x80000000 >> ((i >> 1) & 0x1F)))
			{

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

			}

		}

	}
	
	return 0;
	
}

I have been working on this between class, meals and sleeping.

[axn]

Originally Posted by ShiningArcanine

That does not work, as it sets all of the bits as equal to zero, except the bit being considered. However, your next suggestion will enable me, to use:

arr[i >> 6] ^= ~(0x80000000 >> ((i >> 1) & 0x1F))

It should work. arr[i >> 6] &= ~(0x80000000 >> ((i >> 1) & 0x1F)). This is the standard way of clearing a bit. The ~ will take the NOT of the bit mask and use that to "and" out that bit. No need to check if the bit was already set or not.

Originally Posted by ShiningArcanine

I think that j needs to be used slightly differently because this program does not have memory structures to test even integers.

Already taken care in that code. That's why we start with j = i/2 and increment j by num (j += num) rather than (i += num << 1).

[ShiningArcanine]

Originally Posted by axn1

It should work. arr[i >> 6] &= ~(0x80000000 >> ((i >> 1) & 0x1F)). This is the standard way of clearing a bit. The ~ will take the NOT of the bit mask and use that to "and" out that bit. No need to check if the bit was already set or not.



Already taken care in that code. That's why we start with j = i/2 and increment j by num (j += num) rather than (i += num << 1).

That won't work, as & maps to AND and not XAND, such that if you have 0101 and you want to unset 0100, then ~x = 1010 and ~x AND 0100 = 0000.

If you use XAND, then (~x) XAND 0100 = 0001, but if you were to do it again, you would get what we had originally, as ~x will be 1110, and ~x XAND 0100 = 0101.

(NOT x) XAND y is both a way of setting bits and a way unsetting bits, with the bits you want to either set or unset in x set to 1 in y. That prevents it from being able to replace the conditional expression, as when the sieve goes to 64, the bit position representing 63 will be unset and reset, making the program will report that 63 is prime. However, (NOT x) XAND y == x XOR y, which is x ^ y in C.

I posted the source code. Make your change, compile it in Visual Studio, execute the program, run it for 64 as an upper bound, and watch it report that 63 is prime.

As for your other suggestion, my program only considers odd integers. On the first iteration, j is even because of j = i/2. On the second iteration, j is odd because of (j += num). On the third iteration, j is even because of (j += num). j will oscillate from even to odd over the course of the first inner loop even loop, making the program report false primes and miss primes. Of course, that is ignoring the fact that after the first iteration, the program will sieve for some even numbers, so some loops will consider only even integers and apply their values to bits that represent odd integers.

Edit: I seem to have uncovered a small bug in my program. The inner loop's condition should be i <= len and not i < len. This caused the program to generate incorrect output whenever the upper bound is an odd composite integer. Despite this bug, my 64 upper bound example above still holds. Apparently, I am somewhat biased towards picking either even or prime bounds, so I missed the bug until I accidentally typed 63 instead of 64.

[bsquared]

Originally Posted by ShiningArcanine

That won't work, as & maps to AND and not XAND, such that if you have 0101 and you want to unset 0100, then ~x = 1010 and ~x AND 0100 = 0000.

I'm not familiar with XAND, (exclusive AND?? Wouldn't that always evaluate to 0?) But if you have b = 0101 and want to unset the 0100 then

b &= ~(0100)

works.

[axn]

Originally Posted by ShiningArcanine

That won't work, as & maps to AND and not XAND, such that if you have 0101 and you want to unset 0100, then ~x = 1010 and ~x AND 0100 = 0000.

If you use XAND, then (~x) XAND 0100 = 0001, but if you were to do it again, you would get what we had originally, as ~x will be 1110, and ~x XAND 0100 = 0101.

I think you've confused which pattern is to be NOT-ed.

If you want to unset 0100 of the pattern 0101, then you do 0101 & ~(0100)
i.e. 0101 & 1011 = 0001. I am not sure why you are getting wrong results. Neither XAND nor a test for checking if the bit is set is needed here.
EDIT:- Perhaps, you were checking 63rd bit which corresponds to the prime 127?

Originally Posted by ShiningArcanine

I posted the source code. Make your change, compile it in Visual Studio, execute the program, run it for 64 as an upper bound, and watch it report that 63 is prime.

Don't have access to a C compiler

Originally Posted by ShiningArcanine

As for your other suggestion, my program only considers odd integers. On the first iteration, j is even because of j = i/2. On the second iteration, j is odd because of (j += num). On the third iteration, j is even because of (j += num). j will oscillate from even to odd over the course of the first inner loop even loop, making the program report false primes and miss primes.

The purpose of j is not to report primes; it is used to traverse the bitmap. Clearing the j-th bit of the bitmap is equivalent to clearing the integer 2*j+1. The usage of j is to merely reduce some shifts inside the critical inner loop.

[petrw1]

Originally Posted by wblipp

You could download them and see.

http://www.rsok.com/~jrm/printprimes.html

Also note their commentary on using each bit of an 8-bit byte to respresent the primes in a range of 30 integers. I used this in my program to reduce the size of the disk file. For example even if only one of the 30 integers in a range 3,000,000,001 to 3,000,000,030 is prime you still need 4 bytes to store that number and 32 bytes if all 8 possible primes are prime. With the 8-bit method you always only need 1 byte.