From 28e4f80acfc7c3b5a385356b06edb6b27a04c940 Mon Sep 17 00:00:00 2001 From: laniakea Date: Fri, 10 Apr 2026 00:17:00 +0200 Subject: [PATCH] push --- flake.nix | 30 ++++++++ main.c | 213 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 243 insertions(+) create mode 100644 flake.nix create mode 100644 main.c diff --git a/flake.nix b/flake.nix new file mode 100644 index 0000000..075bd54 --- /dev/null +++ b/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 ]; + }; + }; +} diff --git a/main.c b/main.c new file mode 100644 index 0000000..33181ad --- /dev/null +++ b/main.c @@ -0,0 +1,213 @@ +#include +#include +#include +#include +#include +#include +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; +}