From f7da49b0c0241b05d6bccb7bf68838763a899637 Mon Sep 17 00:00:00 2001
From: "Peter J. Holzer" <hjp@hjp.at>
Date: Fri, 25 Dec 2020 19:45:17 +0100
Subject: [PATCH] Compute primes for 1 second in C

---
 C/Makefile            |  3 +++
 C/eratosthenes_1sec.c | 63 +++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 66 insertions(+)
 create mode 100644 C/Makefile
 create mode 100644 C/eratosthenes_1sec.c

diff --git a/C/Makefile b/C/Makefile
new file mode 100644
index 0000000..8835bd5
--- /dev/null
+++ b/C/Makefile
@@ -0,0 +1,3 @@
+CFLAGS = -O9 -Wall
+
+all: eratosthenes_1sec
diff --git a/C/eratosthenes_1sec.c b/C/eratosthenes_1sec.c
new file mode 100644
index 0000000..1527abf
--- /dev/null
+++ b/C/eratosthenes_1sec.c
@@ -0,0 +1,63 @@
+#include <assert.h>
+#include <inttypes.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+
+typedef uint64_t prime_t;
+
+int main(void) {
+
+    struct timespec ts;
+    double t0, t1;
+
+    clock_gettime(CLOCK_MONOTONIC, &ts);
+    t0 = ts.tv_sec + ts.tv_nsec * 1E-9;
+
+    // Initialize
+    uint8_t *is_prime;
+    prime_t n = 4;
+
+    is_prime = malloc(n);
+    is_prime[0] = 0;
+    is_prime[1] = 0;
+    is_prime[2] = 1;
+    is_prime[3] = 1;
+
+    for (;;) {
+        clock_gettime(CLOCK_MONOTONIC, &ts);
+        t1 = ts.tv_sec + ts.tv_nsec * 1E-9;
+        if (t1 - t0 >= 1.0) break;
+
+        prime_t n1 = n * 2;
+        assert(n1 > n);
+        is_prime = realloc(is_prime, n1);
+        assert(is_prime);
+        memset(is_prime + n, 1, n1 - n);
+
+        for (prime_t i = 2; i < n1; i++) {
+            if (!is_prime[i]) continue;
+            // printf("%ju is prime:", (uintmax_t)i);
+
+            prime_t m = (n + (i - 1)) / i * i;
+            if (m == i) m += i;
+
+            for (; m < n1; m += i) {
+                // printf(" %ju", (uintmax_t)m);
+                is_prime[m] = 0;
+            }
+            // printf("\n");
+        }
+        n = n1;
+    }
+    prime_t n_primes = 0;
+    for (prime_t i = 0; i < n; i++) {
+        if (is_prime[i]) {
+            // printf("%ju\n", (uintmax_t)i);
+            n_primes++;
+        }
+    }
+    printf("%ju primes <= %ju in %g seconds\n", (uintmax_t)n_primes, (uintmax_t)n, t1 - t0);
+
+}