new file: gitignore

new file:   primes/src/sive.c

gitignore

@@ -0,0 +1,6 @@
+*
+!**/*.c
+!**/*.h
+# !**/**/*.h
+# !**/**/*.c
+!**/Makefile

primes/Makefile

@@ -0,0 +1,8 @@
+CFALGSD=-Wall -Wextra -g -ggdb 
+CFALGSF=-Wall -Wextra -Ofast -march=native -mtune=native -funroll-all-loops
+CC = gcc  
+
+fast: src/sive.c 
+	$(CC) $(CFALGSF) src/sive.c -o fastOut
+debug: src/sive.c 
+	$(CC) $(CFALGSD) src/sive.c -o debugOut

primes/src/sive.c

@@ -0,0 +1,189 @@
+#include <assert.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+
+typedef uint64_t WORD;
+typedef uint8_t BOOL;
+const BOOL false = 0;
+const BOOL true = ~false;
+
+// const WORD limit = 542;
+// const WORD limit = 1e2;
+const WORD limit = 1e9;
+const WORD nr = 78498;
+
+// Sive o Atkins hit wheel
+const WORD s[] = {1, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53, 59};
+const WORD wheelLen = 16;
+
+int main() {
+  BOOL *isPrime = calloc(limit, sizeof(WORD));
+  assert(isPrime != NULL);
+
+  /*
+  // Initialize the sieve with enough wheels to include limit:
+for n ← 60 × w + x where w ∈ {0,1,...,limit ÷ 60}, x ∈ s:
+    is_prime(n) ← false
+    */
+  for (WORD w = 0; w < limit / 60; w++) {
+    for (WORD i = 0; i < wheelLen; i++) {
+      WORD n = 60 * w + s[i];
+      if (n >= limit)
+        break;
+      isPrime[n] = false;
+    }
+  }
+
+  /*
+   * // Put in candidate primes:
+//   integers which have an odd number of
+//   representations by certain quadratic forms.
+// Algorithm step 3.1:
+for n ≤ limit, n ← 4x²+y² where x ∈ {1,2,...} and y ∈ {1,3,...} // all x's odd
+y's if n mod 60 ∈ {1,13,17,29,37,41,49,53}: is_prime(n) ← ¬is_prime(n)   //
+toggle state
+    */
+  WORD mask1[] = {1, 13, 17, 29, 37, 41, 49, 53};
+  WORD maskLen1 = 8;
+  for (WORD x = 1; 4 * x * x < limit; x += 1) {
+    for (WORD y = 1; 1; y += 2) {
+      WORD n = 4 * x * x + y * y;
+      if (n >= limit) {
+        break;
+      }
+      WORD mod = n % 60;
+      for (WORD i = 0; i < maskLen1; i++) {
+        if (mod == mask1[i]) {
+          isPrime[n] = ~isPrime[n];
+          // printf("3.1 %lu is %s prime\n", n, (isPrime[n]) ? "" : "not");
+          break;
+        }
+      }
+    }
+  }
+
+  /*
+  // Algorithm step 3.2:
+for n ≤ limit, n ← 3x²+y² where x ∈ {1,3,...} and y ∈ {2,4,...} // only odd x's
+    if n mod 60 ∈ {7,19,31,43}:                                 // and even y's
+        is_prime(n) ← ¬is_prime(n)   // toggle state
+    */
+  WORD mask2[] = {7, 19, 31, 43};
+  WORD maskLen2 = 4;
+  for (WORD x = 1; x * x * 3 < limit; x += 2) {
+    for (WORD y = 2; y * y < limit; y += 2) {
+
+      WORD n = 3 * x * x + y * y;
+      if (n >= limit) {
+        break;
+      }
+      WORD mod = n % 60;
+      for (WORD i = 0; i < maskLen2; i++) {
+        if (mod == mask2[i]) {
+          isPrime[n] = ~isPrime[n];
+          // printf("3.2 %lu is %s prime\n", n, (isPrime[n]) ? "" : "not");
+          break;
+        }
+      }
+    }
+  }
+
+  /*
+  // Algorithm step 3.3:
+  for n ≤ limit, n ← 3x²-y² where x ∈ {2,3,...} and y ∈ {x-1,x-3,...,1} //all
+  even/odd if n mod 60 ∈ {11,23,47,59}:                                   //
+  odd/even combos is_prime(n) ← ¬is_prime(n)   // toggle state
+  */
+  WORD mask3[] = {11, 23, 47, 59};
+  WORD maskLen3 = 4;
+  for (WORD x = 2; (3 * x * x) - ((x - 1) * (x - 1)) < limit; x += 1) {
+    WORD y = x - 1;
+    while (y >= 1) {
+      // printf("x: %lu, y:%lu\n", x, y);
+      // sleep(1);
+
+      WORD n = 3 * x * x - y * y;
+
+      if (n >= limit) {
+        // breakOuter = true;
+        if (y < 2) {
+          break;
+        }
+        y -= 2;
+        continue;
+      }
+      WORD mod = n % 60;
+      for (WORD i = 0; i < maskLen3; i++) {
+        if (mod == mask3[i]) {
+          isPrime[n] = ~isPrime[n];
+          // printf("3.3 %lu is %s prime\n", n, (isPrime[n]) ? "" : "not");
+          break;
+        }
+      }
+      y = (y > 2) ? y - 2 : 0;
+    }
+  }
+
+  /*
+   *
+// Eliminate composites by sieving, only for those occurrences on the wheel:
+for n² ≤ limit, n ← 60 × w + x where w ∈ {0,1,...}, x ∈ s, n ≥ 7:
+    if is_prime(n):
+        // n is prime, omit multiples of its square; this is sufficient
+        // because square-free composites can't get on this list
+        for c ≤ limit, c ← n² × (60 × w + x) where w ∈ {0,1,...}, x ∈ s:
+            is_prime(c) ← false
+
+    */
+
+  for (WORD w = 0; w / 60 < limit; w++) {
+    WORD n = ~0;
+    for (WORD i = 0; i < wheelLen; i++) {
+      n = 60 * w + s[i];
+      if (n * n >= limit) {
+        break;
+      }
+      // or break;?
+
+      if (n < limit && isPrime[n]) {
+        for (WORD ww = 0; n * n * (60 * ww) < limit; ww++) {
+          WORD c;
+          for (WORD ii = 0; ii < wheelLen; ii++) {
+            c = n * n * (60 * ww + s[ii]);
+            if (c < limit) {
+              isPrime[c] = false;
+              // printf("elim: %lu is %s prime\n", c, (isPrime[c]) ? "" :
+              // "not");
+            }
+          }
+        }
+      }
+    }
+  }
+
+  // isPrime[2] = true;
+  // isPrime[3] = true;
+  // isPrime[5] = true;
+  printf("2 3 5");
+  WORD sum = 3;
+  for (WORD w = 0; w / 60 < limit; w++) {
+    WORD n = ~0;
+    for (WORD i = 0; i < wheelLen; i++) {
+      n = 60 * w + s[i];
+      if (n < 7)
+        continue;
+      if (n >= limit)
+        break;
+      if (isPrime[n]) {
+        // printf(" %lu ", n);
+        sum += 1;
+      }
+    }
+  }
+  printf("\n");
+  printf("%lu\n", sum);
+
+  return 0;
+}