push

flake.nix

@@ -0,0 +1,30 @@
+{
+  description = "Lucy prime counting sieve";
+  inputs.nixpkgs.url = "github:NixOS/nixpkgs/nixos-unstable";
+
+  outputs = { self, nixpkgs }: let
+    system = "x86_64-linux";
+    pkgs = nixpkgs.legacyPackages.${system};
+  in {
+    packages.${system}.default = pkgs.stdenv.mkDerivation {
+      pname = "lucy";
+      version = "1.0";
+      src = ./.;
+
+      buildInputs = [ pkgs.glibc ];
+
+      buildPhase = ''
+        gcc -O2 -o lucy main.c -lm
+      '';
+
+      installPhase = ''
+        mkdir -p $out/bin
+        cp lucy $out/bin/lucy
+      '';
+    };
+
+    devShells.${system}.default = pkgs.mkShell {
+      packages = [ pkgs.gcc pkgs.gdb pkgs.valgrind ];
+    };
+  };
+}

main.c

@@ -0,0 +1,213 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <string.h>
+#include <math.h>
+#include <stdbool.h>
+typedef struct {
+    int64_t  x;
+    int64_t  isqrt;   
+    int64_t *arr;     
+    int      len;     
+} FIArray;
+
+static int64_t isqrt64(int64_t n) {
+    int64_t r = (int64_t)sqrt((double)n);
+    while (r * r > n) r--;
+    while ((r+1)*(r+1) <= n) r++;
+    return r;
+}
+
+static int fi_len(int64_t x, int64_t sq) {
+
+    if (sq * sq == x) return (int)(2*sq - 1);
+    return (int)(2*sq);
+}
+
+static FIArray new_fi_array(int64_t x) {
+    FIArray fa;
+    fa.x     = x;
+    fa.isqrt = isqrt64(x);
+    fa.len   = fi_len(x, fa.isqrt);
+    fa.arr   = (int64_t*)calloc(fa.len, sizeof(int64_t));
+    return fa;
+}
+
+static void free_fi_array(FIArray *fa) {
+    free(fa->arr);
+    fa->arr = NULL;
+}
+
+static inline int64_t fi_get(const FIArray *fa, int64_t v) {
+    if (v <= fa->isqrt)
+        return fa->arr[v - 1];
+    else
+        return fa->arr[fa->len - (int)(fa->x / v)];
+}
+
+static inline void fi_set(FIArray *fa, int64_t v, int64_t val) {
+    if (v <= fa->isqrt)
+        fa->arr[v - 1] = val;
+    else
+        fa->arr[fa->len - (int)(fa->x / v)] = val;
+}
+
+static inline int64_t fi_tail(const FIArray *fa, int i) {
+    return fa->arr[fa->len - i];
+}
+static inline void fi_tail_sub(FIArray *fa, int i, int64_t delta) {
+    fa->arr[fa->len - i] -= delta;
+}
+
+static int64_t *keys_fi(int64_t x, int64_t sq, int len) {
+    int64_t *V = (int64_t*)malloc(len * sizeof(int64_t));
+    int idx = 0;
+    int64_t small_max = (sq*sq == x) ? sq-1 : sq;
+    for (int64_t v = 1; v <= small_max; v++)
+        V[idx++] = v;
+    int64_t start_i = (sq*sq == x) ? sq : sq;
+    for (int64_t i = start_i; i >= 1; i--)
+        V[idx++] = x / i;
+    return V;
+}
+
+typedef struct {
+    int64_t *tree;
+    int      n;
+} Fenwick;
+
+static Fenwick new_fenwick(int n, int64_t init_val) {
+    Fenwick fw;
+    fw.n    = n;
+    fw.tree = (int64_t*)calloc(n + 1, sizeof(int64_t));
+    if (init_val != 0) {
+        for (int i = 1; i <= n; i++) {
+            fw.tree[i] += init_val;
+            int j = i + (i & -i);
+            if (j <= n) fw.tree[j] += fw.tree[i];
+        }
+    }
+    return fw;
+}
+
+static void free_fenwick(Fenwick *fw) {
+    free(fw->tree);
+    fw->tree = NULL;
+}
+
+static int64_t fen_sum(const Fenwick *fw, int i) {
+    int64_t s = 0;
+    for (; i > 0; i -= i & -i)
+        s += fw->tree[i];
+    return s;
+}
+
+static void fen_add(Fenwick *fw, int i, int64_t delta) {
+    for (; i <= fw->n; i += i & -i)
+        fw->tree[i] += delta;
+}
+
+static void fen_zero(Fenwick *fw, int i) {
+    fen_add(fw, i, -1);  
+}
+
+
+static FIArray lucy(int64_t x) {
+    FIArray S  = new_fi_array(x);
+    int64_t sq = S.isqrt;
+    int64_t *V = keys_fi(x, sq, S.len);
+
+    for (int i = 0; i < S.len; i++)
+        S.arr[i] = V[i] - 1;
+
+    for (int64_t p = 2; p <= sq; p++) {
+        if (S.arr[p-1] == S.arr[p-2]) continue; 
+        int64_t sp = S.arr[p-2]; 
+        for (int i = S.len - 1; i >= 0; i--) {
+            int64_t v = V[i];
+            if (v < p * p) break;
+            S.arr[i] -= fi_get(&S, v / p) - sp;
+        }
+    }
+
+    free(V);
+    return S;
+}
+
+static FIArray lucy_fenwick(int64_t x) {
+    FIArray S  = new_fi_array(x);
+    int64_t sq = S.isqrt;
+    int64_t y;
+    if (x == 1) {
+        y = 1;
+    } else {
+        double xf = (double)x;
+        y = (int64_t)round(0.35 * pow(xf, 2.0/3.0) / pow(log(xf), 2.0/3.0));
+        int64_t ymax = (int64_t)4e9;
+        if (y > ymax) y = ymax;
+        if (y < sq + 1) y = sq + 1;
+    }
+
+    bool    *sieve_raw = (bool*)calloc(y + 1, sizeof(bool));
+    Fenwick  sieve     = new_fenwick((int)y, 1); 
+    sieve_raw[0] = true; fen_add(&sieve, 1, -1); 
+    sieve_raw[1] = true; fen_add(&sieve, 2, -1); 
+
+    fen_add(&sieve, 1, -1); 
+
+    free_fenwick(&sieve);
+    free(sieve_raw);
+
+    sieve_raw = (bool*)calloc(y + 1, sizeof(bool));
+    sieve     = new_fenwick((int)y, 1);
+    sieve_raw[1] = true;
+    fen_add(&sieve, 1, -1);
+
+    int64_t *V = keys_fi(x, sq, S.len);
+    for (int i = 0; i < S.len; i++)
+        S.arr[i] = V[i] - 1;
+
+    #define S0(v) ((v) <= y ? fen_sum(&sieve, (int)(v)) : fi_get(&S, (v)))
+
+    for (int64_t p = 2; p <= sq; p++) {
+        if (sieve_raw[p]) continue; 
+        int64_t sp = fen_sum(&sieve, (int)(p - 1));
+        int64_t lim = x / y;
+        if (x / (p * p) < lim) lim = x / (p * p);
+
+        for (int64_t i = 1; i <= lim; i++) {
+            fi_tail_sub(&S, (int)i, S0(x / (i * p)) - sp);
+        }
+
+        for (int64_t j = p * p; j <= y; j += p) {
+            if (!sieve_raw[j]) {
+                sieve_raw[j] = true;
+                fen_add(&sieve, (int)j, -1);
+            }
+        }
+    }
+
+    for (int i = 0; i < S.len; i++) {
+        if (V[i] > y) break;
+        S.arr[i] = fen_sum(&sieve, (int)V[i]);
+    }
+
+    #undef S0
+    free(V);
+    free(sieve_raw);
+    free_fenwick(&sieve);
+    return S;
+}
+
+int main(void) {
+    int64_t x = (int64_t)1e14;  /* 10^14 */
+
+    printf("Computing pi(%lld) ...\n", (long long)x);
+    FIArray S = lucy_fenwick(x);
+
+    int64_t result = fi_get(&S, x);
+    printf("pi(10^14) = %lld\n", (long long)result);
+
+    free_fi_array(&S);
+    return 0;
+}