Project Euler problem 10

[Pyrofan1]

I'm trying to solve Project Euler problem 10 which is to figure out the sum of all primes under 2 million.
Here's my code

Code: Select all

#include <stdio.h>

int isprime(unsigned int num)
{
	if(num<10)
	{
		if(num==2) return 1;
		if(num==3) return 1;
		if(num==5) return 1;
		if(num==7) return 1;
	}
	else
	if(num%2==0||num%3==0||num%5==0||num%7==0) return 0;
	else return 1;
}

int main()
{
	unsigned int count=0;
	unsigned int ans=0;

	for(count=0;count!=2000000;count++)
	{
		if(isprime(count)==1)
		{
			//printf("%d\n",count);
			ans+=count;
		}
	}

	printf("%d\n",ans);
	return 0;
}

Problem #1: it doesn't compute the right answer
Problem #2: the answer it gives changes depending on if the printf("%d\n",count); line is commented or not

problem #2 is concerning me the most because when the line is commented it gives the answer as 1876609528, but when it's not commented it gives 1876609501 and neither of these answers are correct.

[Shadyjames]

Numbers that aren't divisible by 2, 3, 5 or 7 aren't necessarily prime. Take 143 for example, which is 11*13. Don't know how to help you much more than that without just giving away the answer.
No idea about the printf thing, i don't do C.

[Troy Martin]

http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes

[eddyb]

Code: Select all

int is_not_prime_table[2000000];
int i, i2;
int curr;
memset(&is_not_prime_table, 0, 8000000);//2 millions * sizeof(int)
for(i=1;i<2000000;i++) {
       if(!is_not_prime_table[i]) { //if is a prime
             curr+=i;
             for(i2=1;i2<2000000/i;i++) { //for each multiple of i
                     is_not_prime_table[i*i2]=1; //mark it as not a prime
             }
       }
}
printf("%d", curr);

Hope this would help you (I'm good at maths ).

[Pyrofan1]

I fixed my code and came up with this:

Code: Select all

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

int isprime(uint64_t num)
{
	unsigned int count;
	unsigned int max=(unsigned int)sqrt(num)+1;
	
	if(num<10)
	{
		if(num==2||num==3||num==5||num==7) return 1;
		else return 0;
	}
	else
	{
		for(count=2;count<max;count++)
		{
			if(num%count==0) return 0;
		}
	}
	return 1;
}

int main()
{
	uint64_t count=0;
	uint64_t ans=0;

	for(count=2;count<2000000;count++)
	{
		if(isprime(count)==1)
		{
			ans+=count;
		}
	}

	printf("%u\n",ans);
	return 0;
}

my code gives the answer as 1179908154
Shiner's code gives the answer as 4271563928
And the website says that neither is right.

I'm guessing that Shiner's code is over flowing the int and trying it using uint64_t numbers gives the answer 4140393044 which isn't right either

[JohnnyTheDon]

Your isprime() function doesn't work. It will only find some primes. The easiest way to calculate the primes is the Sieve of Eratosthenes that Troy suggested.

Here is some code:

Code: Select all

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

#define HIGHEST 2000000

char Prime[HIGHEST];

int main(int argc, char** argv)
{
	/* Set entires array to true (prime) */
	memset(Prime, 1, sizeof(char)*HIGHEST);
	/* except for 1 */
	Prime[0] = 0;
	
	char done = 0;
	unsigned int strikeOut = 2;
	
	while(!done)
	{
		unsigned int i;
		for(i = 2; i < HIGHEST; i++)
		{
			if((i+1)%strikeOut == 0 && (i+1) != strikeOut)
				Prime[i] = 0;
		}
		for(i = strikeOut; i < HIGHEST; i++)
		{
			if(Prime[i])
			{
				strikeOut = i+1;
				break;
			}
			else if((double)(i+1) > sqrt(HIGHEST))
			{
				done = 1;
				break;
			}
		}
	}
	
	int i;
	unsigned long long sum = 0;
	for(i = 1; i < HIGHEST; i++)
	{
		if(Prime[i])
		{
			sum += (unsigned long long)(i+1);
		}
	}
	
	printf("Sum: %llu\n", sum);
}

Answer: 142,913,828,922

[Brendan]

Hi,

Here is some code:

Hmm - with "gcc -O3" that took 4.173 seconds on my computer...

I optimized:

Code: Select all

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

#define HIGHEST 2000000

char Prime[HIGHEST / 2];

int main(int argc, char** argv)
{
	unsigned int i, j;
	unsigned long long sum = 2;
	unsigned int total = 0;

	/* Set entire array to true (prime) */
	memset(Prime, 1, sizeof(char)*HIGHEST / 2);
	/* except for 1 */
	Prime[0] = 0;

	for(i = 3; i < HIGHEST; i += 2) {
		if(Prime[i / 2] == 1) {
			sum += (unsigned long long)i;
			total++;
		}
		for(j = (i+i+i)/2; j < HIGHEST/2; j += i) {
			Prime[j] = 0;
		}
	}
	printf("Sum: %llu (%d prime numbers)\n", sum, total);
}

That only takes 0.021 seconds on my computer....

If I wasn't lazy, I'd use a bitfield (and the BT and BTR instructions) to improve cache usage (so it only works on 125000 bytes instead of 1000000 bytes), but...


Cheers,

Brendan

[eddyb]

Answer: 142,913,828,922

Well, the answer is good. but what i've wrong?! my code was marking those numbers that are not primes(his multiples) and add it when finds a prime number.

[kasper]

Hi all,

Just out of curiosity. I wonder how the printf of count, in the original code, can affect the answer.

Original code repeated:

[...]

Code: Select all

#include <stdio.h>

int isprime(unsigned int num)
{
	if(num<10)
	{
		if(num==2) return 1;
		if(num==3) return 1;
		if(num==5) return 1;
		if(num==7) return 1;
	}
	else
	if(num%2==0||num%3==0||num%5==0||num%7==0) return 0;
	else return 1;
}

int main()
{
	unsigned int count=0;
	unsigned int ans=0;

	for(count=0;count!=2000000;count++)
	{
		if(isprime(count)==1)
		{
			//printf("%d\n",count);
			ans+=count;
		}
	}

	printf("%d\n",ans);
	return 0;
}

[...]
Problem #2: the answer it gives changes depending on if the printf("%d\n",count); line is commented or not

problem #2 is concerning me the most because when the line is commented it gives the answer as 1876609528, but when it's not commented it gives 1876609501 and neither of these answers are correct.

Any ideas what happens?

Kasper

[Brendan]

Hi,

Well, the answer is good. but what i've wrong?! my code was marking those numbers that are not primes(his multiples) and add it when finds a prime number.

Your method should be even faster, as there's no need to mark non-primes. Your code didn't compile though, and had a few other bugs - I changed it to this:

Code: Select all

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

char is_not_prime_table[2000000];
unsigned int i, i2;
unsigned long long curr = 0;

int main(int argc, char** argv) {
	memset(is_not_prime_table, 0, 2000000);
	is_not_prime_table[0] = 1;
	is_not_prime_table[1] = 1;

	for(i = 2; i < 2000000; i++) {
      	if(is_not_prime_table[i] == 0) {
			curr += (unsigned long long)i;
			for(i2 = i + i; i2 < 2000000; i2 += i) {
				is_not_prime_table[i2] = 1;
			}
		}
	}
	printf("%llu\n", curr);
}

With "gcc -O3" that takes 0.021 seconds on my computer - the same as my optimized version.

Of course I optimized your version by skipping every multiple of 2, and got this:

Code: Select all

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

char is_not_prime_table[1000000];
unsigned int i, i2;
unsigned long long curr = 2;

int main(int argc, char** argv) {
	memset(is_not_prime_table, 0, 1000000);

	for(i = 1; i < 1000000; i++) {
      	if(is_not_prime_table[i] == 0) {
			curr += (unsigned long long)i + (unsigned long long)i + 1;
			for(i2 = i + i + 1 + i; i2 < 1000000; i2 += i + i + 1) {
				is_not_prime_table[i2] = 1;
			}
		}
	}
	printf("%llu\n", curr);
}

That takes half as long (0.010 seconds on my computer)...

Then I did a version in assembly with BT and BTS instructions. It didn't help (took 0.12 seconds on my computer), probably because my CPU has 4 MiB L2 caches and the previous versions don't end up thrashing the cache anyway. That annoyed me, so I sunk my teeth in...

In the end I came up with this:

Code: Select all

section .data
	align 8
primes:	times (1000000) db 0
section .text

getPrimes:
	push rbx
	push rcx
	push rdx

	xor ebx,ebx
	mov eax,2
	jmp .prime1A

.nextGroupA:
	cmp qword [primes+rbx],0x0101010101010101
	je .doneGroupA

	cmp byte [primes+rbx],0
	jne .prime1A
	lea rdx,[rbx*2 + 0*2 + 1]
	lea ecx,[primes+ebx*3 + 0*3 + 1]
	add rax,rdx
	cmp ecx,primes+1000000
	jae .prime1A
.nextMark0:
	mov byte [ecx],1
	add ecx,edx
	cmp ecx,primes+1000000
	jb .nextMark0

.prime1A:
	cmp byte [primes+rbx+1],0
	jne .prime2A
	lea rdx,[rbx*2 + 1*2 + 1]
	lea ecx,[primes+ebx*3 + 1*3 + 1]
	add rax,rdx
	cmp ecx,primes+1000000
	jae .prime2A
.nextMark1:
	mov byte [ecx],1
	add ecx,edx
	cmp ecx,primes+1000000
	jb .nextMark1

.prime2A:
	cmp byte [primes+rbx+2],0
	jne .prime3A
	lea rdx,[rbx*2 + 2*2 + 1]
	lea ecx,[primes+ebx*3 + 2*3 + 1]
	add rax,rdx
	cmp ecx,primes+1000000
	jae .prime3A
.nextMark2:
	mov byte [ecx],1
	add ecx,edx
	cmp ecx,primes+1000000
	jb .nextMark2

.prime3A:
	cmp byte [primes+rbx+3],0
	jne .prime4A
	lea rdx,[rbx*2 + 3*2 + 1]
	lea ecx,[primes+ebx*3 + 3*3 + 1]
	add rax,rdx
	cmp ecx,primes+1000000
	jae .prime4A
.nextMark3:
	mov byte [ecx],1
	add ecx,edx
	cmp ecx,primes+1000000
	jb .nextMark3

.prime4A:
	cmp byte [primes+rbx+4],0
	jne .prime5A
	lea rdx,[rbx*2 + 4*2 + 1]
	lea ecx,[primes+ebx*3 + 4*3 + 1]
	add rax,rdx
	cmp ecx,primes+1000000
	jae .prime5A
.nextMark4:
	mov byte [ecx],1
	add ecx,edx
	cmp ecx,primes+1000000
	jb .nextMark4

.prime5A:
	cmp byte [primes+rbx+5],0
	jne .prime6A
	lea rdx,[rbx*2 + 5*2 + 1]
	lea ecx,[primes+ebx*3 + 5*3 + 1]
	add rax,rdx
	cmp ecx,primes+1000000
	jae .prime6A
.nextMark5:
	mov byte [ecx],1
	add ecx,edx
	cmp ecx,primes+1000000
	jb .nextMark5

.prime6A:
	cmp byte [primes+rbx+6],0
	jne .prime7A
	lea rdx,[rbx*2 + 6*2 + 1]
	lea ecx,[primes+ebx*3 + 6*3 + 1]
	add rax,rdx
	cmp ecx,primes+1000000
	jae .prime7A
.nextMark6:
	mov byte [ecx],1
	add ecx,edx
	cmp ecx,primes+1000000
	jb .nextMark6

.prime7A:
	cmp byte [primes+rbx+7],0
	jne .doneGroupA
	lea rdx,[rbx*2 + 7*2 + 1]
	lea ecx,[primes+ebx*3 + 7*3 + 1]
	add rax,rdx
	cmp ecx,primes+1000000
	jae .doneGroupA
.nextMark7:
	mov byte [ecx],1
	add ecx,edx
	cmp ecx,primes+1000000
	jb .nextMark7

.doneGroupA:
	add ebx,8
	cmp ebx,(1000000-1)/3+7
	jb .nextGroupA

	xor edx,edx
.nextGroupB:
	cmp qword [primes+rbx],0x0101010101010101
	je .doneGroupB

	cmp byte [primes+rbx],0
	jne .prime1B
	lea rax,[rax+rbx*2 + 0*2 + 1]
.prime1B:
	cmp byte [primes+rbx+1],0
	jne .prime2B
	lea rdx,[rdx+rbx*2 + 1*2 + 1]
.prime2B:
	cmp byte [primes+rbx+2],0
	jne .prime3B
	lea rax,[rax+rbx*2 + 2*2 + 1]
.prime3B:
	cmp byte [primes+rbx+3],0
	jne .prime4B
	lea rdx,[rdx+rbx*2 + 3*2 + 1]
.prime4B:
	cmp byte [primes+rbx+4],0
	jne .prime5B
	lea rax,[rax+rbx*2 + 4*2 + 1]
.prime5B:
	cmp byte [primes+rbx+5],0
	jne .prime6B
	lea rdx,[rdx+rbx*2 + 5*2 + 1]
.prime6B:
	cmp byte [primes+rbx+6],0
	jne .prime7B
	lea rax,[rax+rbx*2 + 6*2 + 1]
.prime7B:
	cmp byte [primes+rbx+7],0
	jne .doneGroupB
	lea rdx,[rdx+rbx*2 + 7*2 + 1]

.doneGroupB:
	add ebx,8
	cmp ebx,1000000
	jb .nextGroupB

	add rax,rdx

	pop rdx
	pop rcx
	pop rbx
	ret

That takes 0.008 seconds on my computer - roughly 20% faster than C...


Cheers,

Brendan

[JohnnyTheDon]

Just out of curiosity. I wonder how the printf of count, in the original code, can affect the answer.

The reason the printf affected the answer MIGHT be that you are using the format string for a signed integer ( %d ) with an unsigned integer. I don't know why that would affect it, but that is the only reason I can think of. Try switching it to %u

Hmm - with "gcc -O3" that took 4.173 seconds on my computer...

Yeah my code definitely sucks. I was just trying to throw something out that got the right answer.

[Rook]

If you are using MinGW, try "%I64u" instead of "%llu" in printf. This gives me the right answer.

[stephenj]

If you want to generate primes quickly, use the Sieve of Atkin.

There is a 32-bit optimized version at the bottom on the link.

[Rook]

For completeness, here is a fairly fast IsPrime function in C, written following the recommendations on Project Euler:

Code: Select all

// The following code is public domain licenced
#include <math.h>
#include <stdint.h>

uint8_t IsPrime( uint32_t n )
{
    if( n == 2 || n == 3 || n == 5 || n == 7 ) return 1;
    if( n <= 1 || n % 2 == 0 || n % 3 == 0 ) return 0;

    // from here we know that n >= 11 and is odd and not a multiple of 3

    // note: to the right and to the left of every number divisible by 6,
    // i.e. 6k-1 on the left, 6k+1 on the right, is an odd number which
    // is not a multiple of 3

    // if a number is not divisible by some k <= r then it is prime, because
    // no number can have more than one prime factor > floor( its sqrt )

    uint32_t r = (uint32_t)sqrtl( (long double)n );

    for( uint32_t k = 5; k <= r; k += 6 )
    {
        if( n % k == 0 ) return 0;
        if( n % (k+2) == 0 ) return 0;
    }

    return 1;
}

I used it to output the first milion primes into a text file in a matter of seconds.

[Shadyjames]

http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes

Well, i wasn't going to give away the answer unless explicitly asked to do so.

Bit late now though =P