new file: primeCalc/.gitignore

new file:   primeCalc/Makefile
	new file:   primeCalc/README.md
	new file:   primeCalc/src/prime_ex.c
	new file:   primeCalc/src/sieve.c
	new file:   primeCalc/src/sieve.h

primeCalc/.gitignore

@@ -0,0 +1,3 @@
+prime.txt
+fastOut
+debugOut

primeCalc/Makefile

@@ -0,0 +1,18 @@
+CFLAGSD=-g -ggdb -Wall -Wextra -Werror 
+CFLAGSF=-Wall -Wextra -Ofast -march=native -mtune=native -funroll-all-loops
+FILES=src/sieve.c src/sieve.h src/prime_ex.c
+
+all: fast
+
+debug: $(FILES)
+	gcc $(CFLAGSD) $(FILES) -o debugOut 
+
+fast: $(FILES)
+	gcc $(CFLAGSF) $(FILES) -o fastOut 
+
+clean: 
+	rm -rf *.o debugOut fastOut prime.txt
+
+
+
+

primeCalc/README.md

@@ -0,0 +1,10 @@
+# Prime Sieve
+
+Prime sieve using bitfields
+
+'make' for fast benchmark output
+
+'make debug' for debugOut
+
+If program does not finish in appropriate time adjust n_per_clock in src/sieve.c 
+

primeCalc/src/prime_ex.c

@@ -0,0 +1,22 @@
+#include "sieve.h"
+#include <stdio.h>
+#include <time.h>
+
+int main(int argc, char **argv) {
+  FILE *ptr_file;
+  ptr_file = fopen("prime.txt", "w");
+
+  clock_t t_0 = clock();
+  //// prime algorithm START ////
+  sieve(argc, argv, ptr_file);
+
+  //// prime algorithm END ////
+  clock_t t_E = clock();
+
+  fclose(ptr_file);
+
+  float t_delta = (float)(t_E - t_0) / CLOCKS_PER_SEC;
+  printf("\nTime in Seconds: %.3f\n", t_delta);
+
+  return 0;
+}

primeCalc/src/sieve.c

@@ -0,0 +1,109 @@
+/*Author: Simon Schurti
+ *License: MIT
+ *
+ *
+ * Prime sieve using bitfields & masks
+ */
+
+#include "sieve.h"
+
+// primes up to n per clock
+// may need adjusting on different platform to finish in appropiate time
+const WORD n_per_clock = 130;
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+
+WORD upperLimit = 4e9;
+// WORD upperLimit = 542;
+
+// fast square root approximation function for unsigned ints
+//  uses Herons Method
+//  https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Heron's_method
+WORD root(WORD n) {
+
+  WORD x = 0x1 << (MSB(n) - 2); // initial estimate 2 orders smaller than n
+  WORD prev = 0;
+  while (x != prev) {
+    prev = x;
+    x = (x + (n / x)) / 2;
+  }
+  return x;
+}
+
+// find next set bit > prime
+// this will be next prime as all numbers up to next prime should be eliminated
+WORD findNextPrime(WORD *sieve, WORD prime) {
+  WORD next = prime + 2;
+  /*printf("word: %lb\n", sieve[WordIdx(next)]);*/
+  while (IsSet(sieve[WordIdx(next)], next) == 0) {
+    // printf("%lu is not prime\n", next);
+    next += 2;
+    if (next >= upperLimit)
+      return 0;
+  }
+  return next;
+}
+
+int sieve(int argc, char **argv, FILE *fp) {
+  int targetTime = 15;
+  if (argc > 1) {
+    targetTime = atoi(argv[1]);
+  }
+  upperLimit = (targetTime * CLOCKS_PER_SEC) * n_per_clock;
+
+  // upperLimit = 1e6;
+
+  printf("calculating0 up to : %lu\n", upperLimit);
+
+  // 1 bit per prime candidate / 2 beacause all even can be ignored
+  WORD buffSize = upperLimit / 2 / BitsPerWord;
+
+  // main sieve buffer
+  WORD *sieve = (WORD *)calloc(buffSize, BitsPerWord);
+  if (sieve == NULL) {
+    perror("calloc: ");
+    exit(1);
+  }
+
+  // initialize all as prime
+
+  WORD i = buffSize + 1;
+  while (i--) {
+    sieve[i] |= ones;
+  }
+
+  WORD idx;
+  WORD nextPrime = 3; // 2 can be assumed 3=2+1
+  WORD MaxPrime = nextPrime;
+  WORD primesFound = 1; // 2 is prime
+
+  //// findNextPrime returns 0 if none found in limit
+  while (nextPrime > 0 && nextPrime < upperLimit) {
+    primesFound++;
+    WORD maxFactor = upperLimit / nextPrime;
+    for (WORD factor = nextPrime; factor <= maxFactor; factor += 2) {
+      idx = factor * nextPrime;
+      if (idx >= upperLimit)
+        break;
+      /*printf("eliminating %lu because %lu * %lu = %lu\n", idx, factor,*/
+      /*       nextPrime, factor * nextPrime);*/
+      // printf("with mask: %llb\n", BitIdxMask(idx));
+
+      sieve[WordIdx(idx)] &=
+          sieve[WordIdx(idx)] ^
+          (BitIdxMask(idx)); // set multiple of next prime to 0
+
+      // printf("new bu[%lu]: %lb\n", WordIdx(idx), sieve[WordIdx(idx)]);
+    }
+    // printf("-----> %lu  %llb\n", nextPrime, BitIdxMask(nextPrime));
+    MaxPrime = nextPrime;
+    fprintf(fp, "%lu\n", nextPrime); // write prime number to textfile
+    nextPrime = findNextPrime(sieve, nextPrime);
+  }
+  printf("found %lu primes\nmax prime: %lu\n", primesFound, MaxPrime);
+  free(sieve);
+
+  return 0;
+}

primeCalc/src/sieve.h

@@ -0,0 +1,28 @@
+#ifndef SIEVE
+#define SIEVE
+
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#define WORD uint64_t
+#define BYTE uint8_t
+#define BitsPerWord 64
+#define ones 0xFFFFFFFFFFFFFFFFuLL
+
+// find index of word that contains nth bit in total buffer /2 because no evens
+#define WordIdx(n) (((n) / 2) / BitsPerWord)
+// bit offset in word
+#define BitIdx(n) ((((n) / 2) % BitsPerWord))
+
+// mask(3) = 0b10
+#define BitIdxMask(n) 0x1ull << (BitIdx(n))
+
+#define IsSet(w, b) (w & BitIdxMask(b))
+
+// find most significant bit
+//__builtin_clz should work with gcc on all platforms
+#define MSB(n) (BitsPerWord - 1 - __builtin_clz(n))
+int sieve(int argc, char **argv, FILE *fp);
+
+#endif // !SIEVE