From 90d859bd8b79bc7d4cd7228956ae4b5c97b4f60a Mon Sep 17 00:00:00 2001
From: Reiner Herrmann <reiner@reiner-h.de>
Date: Wed, 22 Sep 2010 21:02:15 +0200
Subject: implemented prime sieve in C (_much_ faster); solved projecteuler 087

---
 src/projecteuler/087.c    |  44 +++++++++++++++++++
 src/projecteuler/common.c | 107 ++++++++++++++++++++++++++++++++++++++++++++++
 src/projecteuler/common.h |   7 +++
 3 files changed, 158 insertions(+)
 create mode 100644 src/projecteuler/087.c
 create mode 100644 src/projecteuler/common.c
 create mode 100644 src/projecteuler/common.h

diff --git a/src/projecteuler/087.c b/src/projecteuler/087.c
new file mode 100644
index 0000000..ef12d9f
--- /dev/null
+++ b/src/projecteuler/087.c
@@ -0,0 +1,44 @@
+#include "common.h"
+#include <stdlib.h>
+#include <string.h>
+
+// TODO: currently only works on 64 bit systems
+
+const unsigned long limit = 50000000;
+
+int main(void)
+{
+	unsigned long* p = primes(8000);
+	unsigned long count = p[0];
+	unsigned long i2, i3, i4, counter;
+	char* numbers = malloc(limit*sizeof(char));
+	memset(numbers, 0, limit);
+
+	for(i2=1; i2<=count; i2++)
+		for(i3=1; i3<=count; i3++)
+			for(i4=1; i4<=count; i4++)
+			{
+				unsigned long p4 = p[i4];
+				unsigned long p3 = p[i3];
+				unsigned long p2 = p[i2];
+				unsigned long long result = p4*p4*p4*p4;
+				if(result >= limit) break;
+				result += p3*p3*p3;
+				if(result >= limit) break;
+				result += p2*p2;
+				if(result >= limit) break;
+
+				numbers[result] = 1;
+			}
+
+	counter = 0;
+	for(i2=0; i2<limit; i2++)
+		if(numbers[i2] == 1)
+			counter++;
+
+	printf("%li\n", counter);
+
+	free(p);
+	return 0;
+}
+
diff --git a/src/projecteuler/common.c b/src/projecteuler/common.c
new file mode 100644
index 0000000..15ac8aa
--- /dev/null
+++ b/src/projecteuler/common.c
@@ -0,0 +1,107 @@
+#include <stdlib.h>
+#include <string.h>
+
+extern unsigned long* primes(unsigned long limit)
+{
+	unsigned long i = 2; // start sieving with 2
+	unsigned long pos, x;
+	char* numbers = (char*) malloc(limit*sizeof(char));
+	unsigned long* primes;
+	unsigned long primecount = 0;
+	memset(numbers, 1, limit);
+	numbers[0] = 0;
+
+	while(i*i <= limit)
+	{
+		// sieve current prime
+		x = i*i;
+		while(x <= limit)
+		{
+			numbers[x-1] = 0;
+			x += i;
+		}
+		
+		// search next prime
+		pos = i;
+		while(++pos <= limit)
+		{
+			if(numbers[pos-1] == 1)
+				break;
+		}
+		i = pos;
+	}
+
+	// count primes
+	pos = 0;
+	while(++pos <= limit)
+	{
+		if(numbers[pos-1] == 1)
+			primecount++;
+	}
+
+	// create list of primes
+	primes = (unsigned long*) malloc((primecount+1)*sizeof(unsigned long)); // store count at first element
+	primes[0] = primecount;
+	pos = 0;
+	x = 1;
+	while(++pos <= limit)
+	{
+		if(numbers[pos-1] == 1)
+			primes[x++] = pos;
+	}
+
+	free(numbers);
+	return primes;
+}
+
+extern int isprime(unsigned long number, const unsigned long* primes)
+{
+	unsigned long count = primes[0];
+	unsigned long pos = count/2 + 1;
+	unsigned long dist = pos/2 + 1;
+	unsigned long minpos = 0, maxpos = count;
+
+	if(number > primes[count] || number < 2)
+		return 0;
+
+	// some kind of binary search...
+	while(1)
+	{
+		if(number > primes[pos])
+		{
+			minpos = pos;
+			pos += dist;
+		}
+		else if(number < primes[pos])
+		{
+			maxpos = pos;
+			pos -= dist;
+		}
+		else
+			return 1;
+
+		dist /= 2;
+
+		if(dist == 1)
+		{
+			while(++minpos <= maxpos)
+				if(primes[minpos] == number)
+					return 1;
+			break;
+		}
+	}
+
+	return 0;
+}
+
+/*int main()
+{
+	long* p = primes(10000);
+	//printf("%i\n", isprime(9973, p));
+	int i;
+	for(i=0; i<p[0]; i++)
+		printf("%d: %d\n", p[i+1], isprime(p[i+1], p));
+	free(p);
+	return 0;
+}*/
+
diff --git a/src/projecteuler/common.h b/src/projecteuler/common.h
new file mode 100644
index 0000000..5889bb8
--- /dev/null
+++ b/src/projecteuler/common.h
@@ -0,0 +1,7 @@
+#ifndef COMMON_H
+#define COMMON_H
+
+unsigned long* primes(unsigned long limit);
+int isprime(unsigned long number, const unsigned long* primes);
+
+#endif // COMMON_H
-- 
cgit v1.2.3