Line data Source code
1 : #include <stddef.h>
2 : #include <stdlib.h>
3 : #include <string.h>
4 : #include <inttypes.h>
5 : #include <math.h>
6 : #include <stdio.h>
7 : #include <unistd.h>
8 : #include <pthread.h>
9 :
10 : #include <nut/debug.h>
11 : #include <nut/modular_math.h>
12 : #include <nut/factorization.h>
13 : #include <nut/sieves.h>
14 :
15 6 : uint64_t nut_max_prime_divs(uint64_t max){
16 : //2*3*5*7*11*13*17*19*23*29*31*37*41*43*47
17 6 : if(max < 2ull*3){//the compiler can optimize this to a binary search if it is faster
18 : return 1;
19 6 : }else if(max < 2ull*3*5){
20 : return 2;
21 6 : }else if(max < 2ull*3*5*7){
22 : return 3;
23 6 : }else if(max < 2ull*3*5*7*11){
24 : return 4;
25 4 : }else if(max < 2ull*3*5*7*11*13){
26 : return 5;
27 4 : }else if(max < 2ull*3*5*7*11*13*17){
28 : return 6;
29 4 : }else if(max < 2ull*3*5*7*11*13*17*19){
30 : return 7;
31 0 : }else if(max < 2ull*3*5*7*11*13*17*19*23){
32 : return 8;
33 0 : }else if(max < 2ull*3*5*7*11*13*17*19*23*29){
34 : return 9;
35 0 : }else if(max < 2ull*3*5*7*11*13*17*19*23*29*31){
36 : return 10;
37 0 : }else if(max < 2ull*3*5*7*11*13*17*19*23*29*31*37){
38 : return 11;
39 0 : }else if(max < 2ull*3*5*7*11*13*17*19*23*29*31*37*41){
40 : return 12;
41 0 : }else if(max < 2ull*3*5*7*11*13*17*19*23*29*31*37*41*43){
42 : return 13;
43 0 : }else if(max < 2ull*3*5*7*11*13*17*19*23*29*31*37*41*43*47){
44 0 : return 14;
45 : }
46 : return 15;
47 : }
48 :
49 193 : uint64_t nut_max_primes_le(uint64_t max){
50 193 : return 1.25506*max/log(max);
51 : }
52 :
53 1 : void *nut_sieve_factorizations(uint64_t max, uint64_t *_w){
54 1 : uint64_t w = nut_max_prime_divs(max);
55 1 : size_t pitch = nut_get_factorizations_pitch(w);
56 1 : void *buf = calloc(pitch, max + 1);
57 1 : if(!buf){
58 : return NULL;
59 : }
60 1 : *_w = w;
61 1000000 : for(uint64_t n = 2; n <= max && n; ++n){//TODO: optimize loop
62 999999 : nut_Factors *factors = nut_Pitcharr_get(buf, pitch, n);
63 999999 : if(factors->num_primes){
64 921501 : continue;
65 : }
66 157232 : for(uint64_t p = n, e = 1; p <= max; ++e){
67 3705353 : for(uint64_t m = p; m <= max;){
68 3626619 : nut_Factors *factors = nut_Pitcharr_get(buf, pitch, m);
69 3626619 : uint64_t i = factors->num_primes;
70 3626619 : if(e == 1){
71 2853708 : ++factors->num_primes;
72 : }else{
73 772911 : --i;
74 : }
75 3626619 : factors->factors[i].prime = n;
76 3626619 : factors->factors[i].power = e;
77 3626619 : if(__builtin_add_overflow(m, p, &m)){
78 : break;
79 : }
80 : }
81 78734 : if(__builtin_mul_overflow(p, n, &p)){
82 : break;
83 : }
84 : }
85 : }
86 : return buf;
87 : }
88 :
89 1 : uint64_t nut_get_factorizations_pitch(uint64_t w){
90 1 : static const nut_Factors dummy;
91 1 : return offsetof(nut_Factors, factors) + w*sizeof(dummy.factors[0]);
92 : }
93 :
94 1 : void *nut_sieve_factors(uint64_t max, uint64_t *_w){
95 1 : uint64_t w = nut_max_prime_divs(max);
96 1 : size_t pitch = nut_get_factors_pitch(w);
97 1 : void *buf = calloc(max + 1, pitch);
98 1 : if(!buf){
99 : return NULL;
100 : }
101 1 : *_w = w;
102 1000000 : for(uint64_t n = 2; n <= max && n; ++n){//TODO: optimize loop
103 999999 : nut_u64_Pitcharr *factors = nut_Pitcharr_get(buf, pitch, n);
104 999999 : if(factors->len){
105 921501 : continue;
106 : }
107 2932206 : for(uint64_t m = n; m <= max;){
108 2853708 : nut_u64_Pitcharr *factors = nut_Pitcharr_get(buf, pitch, m);
109 2853708 : factors->elems[factors->len++] = n;
110 2853708 : if(__builtin_add_overflow(m, n, &m)){
111 : break;
112 : }
113 : }
114 : }
115 : return buf;
116 : }
117 :
118 1 : uint8_t *nut_sieve_omega(uint64_t max){
119 1 : uint8_t *buf = calloc(max + 1, sizeof(uint8_t));
120 1 : if(!buf){
121 : return NULL;
122 : }
123 1000 : for(uint64_t n = 2; n <= max && n; ++n){
124 999 : if(buf[n]){
125 831 : continue;
126 : }
127 2294 : for(uint64_t m = n; m <= max;){
128 2126 : buf[m]++;
129 2126 : if(__builtin_add_overflow(m, n, &m)){
130 : break;
131 : }
132 : }
133 : }
134 : return buf;
135 : }
136 :
137 4 : uint64_t *nut_sieve_largest_factors(uint64_t max){
138 4 : uint64_t *buf = calloc(max + 1, sizeof(uint64_t));
139 4 : if(!buf){
140 : return NULL;
141 : }
142 1011219 : for(uint64_t n = 2; n <= max && n; ++n){
143 1011215 : if(buf[n]){
144 931273 : continue;
145 : }
146 2960485 : for(uint64_t m = n; m <= max;){
147 2880543 : buf[m] = n;
148 2880543 : if(__builtin_add_overflow(m, n, &m)){
149 : break;
150 : }
151 : }
152 : }
153 : return buf;
154 : }
155 :
156 1 : uint32_t *nut_sieve_smallest_factors(uint64_t max){
157 1 : uint32_t *buf = calloc(max + 1, sizeof(uint32_t));
158 1 : if(!buf){
159 : return NULL;
160 : }
161 1 : uint64_t rmax = nut_u64_nth_root(max, 2);
162 31 : for(uint64_t n = 2; n <= rmax; ++n){
163 30 : if(buf[n]){
164 19 : continue;
165 : }
166 1422 : for(uint64_t m = n*n; m <= max;){
167 1411 : if(buf[m] == 0 || buf[m] > n){
168 831 : buf[m] = n;
169 : }
170 1411 : if(__builtin_add_overflow(m, n, &m)){
171 : break;
172 : }
173 : }
174 11 : buf[n] = 1; // this means the first store in the above loop is dead
175 : }
176 970 : for(uint64_t n = rmax + 1; n <= max; ++n){
177 969 : if(!buf[n]){
178 157 : buf[n] = 1;
179 : }
180 : }
181 : return buf;
182 : }
183 :
184 2 : uint32_t *nut_sieve_smallest_factors_wheel6(uint64_t max){
185 2 : uint64_t r = max%6;
186 2 : switch(r){
187 0 : case 0: max -= 1; r = 5; break;
188 2 : case 2 ... 4: max -= (r-1); r = 1;
189 : }
190 2 : uint64_t qmax = max/6;
191 2 : uint64_t max_idx = 2*qmax + (r == 5);
192 2 : uint32_t *buf = calloc(max_idx + 1, sizeof(uint32_t));
193 2 : if(!buf){
194 : return NULL;
195 : }
196 2 : uint64_t rmax = nut_u64_nth_root(max, 2);
197 : // 6*qn + 1 <= rmax
198 2 : uint64_t idx = 0;
199 3297 : for(uint64_t qn = 0; qn <= (rmax - 1)/6; ++qn){
200 3295 : idx += 1;
201 3295 : uint64_t p = 6*qn + 1;
202 3295 : if(buf[2*qn] == 0 && qn){ // we need to skip p = 1 and composites
203 : // 6*qk + 1 <= max/p
204 49524161 : for(uint64_t qk = qn; qk <= (max/p - 1)/6; ++qk){
205 : /* (6*qn + rn)*(6*qk + rk)
206 : 6*((6*qn + rn)*qk + qn*rk) + rn*rk
207 : */
208 49523046 : uint64_t idx = 2*(p*qk + qn);
209 49523046 : if(buf[idx] == 0 || buf[idx] > p){
210 22864372 : buf[idx] = p;
211 : }
212 : // 6*qk + 5 <= max/p
213 49523046 : if(qk <= (max/p - 5)/6){
214 49522309 : uint64_t idx = 2*(p*qk + 5*qn) + 1;
215 49522309 : if(buf[idx] == 0 || buf[idx] > p){
216 22879805 : buf[idx] = p;
217 : }
218 : }
219 : }
220 : }
221 3295 : idx += 1;
222 3295 : if(qn){
223 3295 : }
224 3295 : buf[2*qn] = 1;
225 3295 : if(buf[2*qn + 1]){
226 2170 : continue;
227 : }
228 1125 : p += 4;
229 1125 : if(p > rmax){
230 : break;
231 : }
232 58683054 : for(uint64_t qk = qn; qk <= (max/p - 1)/6; ++qk){
233 58681929 : if(qk > qn){
234 58680804 : uint64_t idx = 2*(p*qk + qn) + 1;
235 58680804 : if(buf[idx] == 0 || buf[idx] > p){
236 31618493 : buf[idx] = p;
237 : }
238 : }
239 58681929 : if(qk <= (max/p - 5)/6){
240 58681196 : uint64_t idx = 2*(p*qk + 5*qn + 4);
241 58681196 : if(buf[idx] == 0 || buf[idx] > p){
242 31634897 : buf[idx] = p;
243 : }
244 : }
245 : }
246 1125 : buf[2*qn + 1] = 1;
247 : }
248 129791790 : for(; idx <= max_idx; ++idx){
249 129791788 : if(!buf[idx]){
250 20798569 : buf[idx] = 1;
251 : }
252 : }
253 : return buf;
254 : }
255 :
256 1005215 : void nut_fill_factors_from_largest(nut_Factors *restrict out, uint64_t n, const uint64_t largest_factors[restrict static n + 1]){
257 1005215 : out->num_primes = 0;
258 4653844 : for(uint64_t p = largest_factors[n], k = 1; p;){
259 3648629 : n /= p;
260 3648629 : uint64_t q = largest_factors[n];
261 3648629 : if(p == q){
262 779955 : ++k;
263 : }else{
264 2868674 : nut_Factor_append(out, p, k);
265 2868674 : p = q;
266 2868674 : k = 1;
267 : }
268 : }
269 1005215 : }
270 :
271 1000 : void nut_fill_factors_from_smallest(nut_Factors *restrict out, uint64_t n, const uint32_t smallest_factors[restrict static n + 1]){
272 1000 : out->num_primes = 0;
273 3877 : for(uint64_t p = smallest_factors[n], k = 1; p;){
274 2877 : if(p == 1){
275 168 : p = n;
276 : }
277 2877 : n /= p;
278 2877 : uint64_t q = smallest_factors[n];
279 2877 : if(q == 1){
280 831 : q = n;
281 : }
282 2877 : if(p == q){
283 751 : ++k;
284 : }else{
285 2126 : nut_Factor_append(out, p, k);
286 2126 : p = q;
287 2126 : k = 1;
288 : }
289 : }
290 1000 : }
291 :
292 1000 : void nut_fill_factors_from_smallest_wheel6(nut_Factors *restrict out, uint64_t n, const uint32_t smallest_factors[restrict static n/3 + 1]){
293 1000 : out->num_primes = 0;
294 1000 : uint64_t t = __builtin_ctzll(n);
295 1000 : if(t){
296 500 : nut_Factor_append(out, 2, t);
297 500 : n >>= t;
298 : }
299 1000 : t = 0;
300 1498 : while(n%3 == 0){
301 498 : ++t;
302 498 : n /= 3;
303 : }
304 1000 : if(t){
305 333 : nut_Factor_append(out, 3, t);
306 : }
307 1377 : while(n != 1){
308 1293 : uint64_t p = smallest_factors[n/3];
309 1293 : if(p == 1){
310 916 : nut_Factor_append(out, n, 1);
311 916 : break;
312 : }
313 : t = 0;
314 832 : while(n%p == 0){
315 455 : n /= p;
316 455 : t += 1;
317 : }
318 377 : nut_Factor_append(out, p, t);
319 : }
320 1000 : }
321 :
322 1 : uint64_t nut_get_factors_pitch(uint64_t w){
323 1 : return offsetof(nut_u64_Pitcharr, elems) + w*sizeof(uint64_t);
324 : }
325 :
326 : // Works by multiplying power + 1 for all prime factors
327 4 : uint64_t *nut_sieve_sigma_0(uint64_t max){
328 4 : uint64_t *buf = malloc((max + 1)*sizeof(uint64_t));
329 4 : if(!buf){
330 : return NULL;
331 : }
332 1002223 : for(uint64_t i = 1; i <= max && i; ++i){
333 1002219 : buf[i] = 1;
334 : }
335 1002219 : for(uint64_t n = 2; n <= max && n; ++n){//TODO: optimize loop
336 1002215 : if(buf[n] != 1){
337 923334 : continue;
338 : }
339 2937250 : for(uint64_t m = n; m <= max;){
340 2858369 : uint64_t a = m/n, e = 1;
341 3632935 : while(a%n == 0){
342 774566 : a /= n;
343 774566 : ++e;
344 : }
345 2858369 : buf[m] *= e + 1;
346 2858369 : if(__builtin_add_overflow(m, n, &m)){
347 : break;
348 : }
349 : }
350 : }
351 : return buf;
352 : }
353 :
354 : // Works by multiplying (prime**(power+1)-1)/(prime-1) for all prime factors
355 1 : uint64_t *nut_sieve_sigma_1(uint64_t max){
356 1 : uint64_t *buf = malloc((max + 1)*sizeof(uint64_t));
357 1 : if(!buf){
358 : return NULL;
359 : }
360 1000001 : for(uint64_t i = 1; i <= max && i; ++i){
361 1000000 : buf[i] = 1;
362 : }
363 1000000 : for(uint64_t n = 2; n <= max && n; ++n){//TODO: optimize loop
364 999999 : if(buf[n] != 1){
365 921501 : continue;
366 : }
367 2932206 : for(uint64_t m = n; m <= max;){
368 2853708 : uint64_t a = n*n;
369 3626619 : while(m%a == 0){
370 772911 : a *= n;
371 : }
372 2853708 : buf[m] *= (a - 1)/(n - 1);
373 2853708 : if(__builtin_add_overflow(m, n, &m)){
374 : break;
375 : }
376 : }
377 : }
378 : return buf;
379 : }
380 :
381 : // Works by multiplying (prime**((power+1)*e)-1)/(prime**e-1) for all prime factors, where power is the
382 : // power of each prime and e is the power of divisors to sum
383 1 : uint64_t *nut_sieve_sigma_e(uint64_t max, uint64_t e){
384 1 : uint64_t *buf = malloc((max + 1)*sizeof(uint64_t));
385 1 : if(!buf){
386 : return NULL;
387 : }
388 1000001 : for(uint64_t i = 1; i <= max && i; ++i){
389 1000000 : buf[i] = 1;
390 : }
391 1000000 : for(uint64_t n = 2; n <= max && n; ++n){//TODO: optimize loop
392 999999 : if(buf[n] != 1){
393 921501 : continue;
394 : }
395 2932206 : for(uint64_t m = n; m <= max;){
396 2853708 : uint64_t a = n*n;
397 3626619 : while(m%a == 0){
398 772911 : a *= n;
399 : }
400 2853708 : buf[m] *= (nut_u64_pow(a, e) - 1)/(nut_u64_pow(n, e) - 1);
401 2853708 : if(__builtin_add_overflow(m, n, &m)){
402 : break;
403 : }
404 : }
405 : }
406 : return buf;
407 : }
408 :
409 37 : bool nut_u64_make_factorial_tbl(uint64_t k, uint64_t modulus, uint64_t bits, uint64_t max_denom, uint64_t factorials[static bits + k], uint64_t inv_factorials[static max_denom + 1]){
410 37 : factorials[0] = 1;
411 3947 : for(uint64_t i = 1; i < bits + k; ++i){
412 3910 : factorials[i] = factorials[i - 1]*i%modulus;
413 : }
414 37 : memcpy(inv_factorials, factorials, 2*sizeof(uint64_t));
415 3620 : for(uint64_t i = 2; i <= max_denom; ++i){
416 3583 : int64_t inv;
417 3583 : if(nut_i64_egcd(factorials[i], modulus, &inv, NULL) != 1){
418 0 : return false;
419 3583 : }else if(inv < 0){
420 1827 : inv_factorials[i] = inv + (int64_t)modulus;
421 : }else{
422 1756 : inv_factorials[i] = inv;
423 : }
424 : }
425 : return true;
426 : }
427 :
428 19 : static uint64_t *sieve_dk_mod(uint64_t max, uint64_t k, uint64_t modulus){
429 19 : uint64_t *buf = malloc((max + 1)*sizeof(uint64_t));
430 38 : uint8_t *is_c_buf [[gnu::cleanup(cleanup_free)]] = calloc(max/8 + 1, sizeof(uint8_t));
431 19 : uint64_t bits = 64 - __builtin_clzll(max);
432 19 : uint64_t max_denom = bits > k - 1 ? bits : k - 1;
433 38 : uint64_t *factorials [[gnu::cleanup(cleanup_free)]] = malloc((bits + k)*sizeof(uint64_t));
434 38 : uint64_t *inv_factorials [[gnu::cleanup(cleanup_free)]] = malloc((max_denom + 1)*sizeof(uint64_t));
435 19 : if(!buf || !is_c_buf || !factorials || !inv_factorials || modulus > INT64_MAX){
436 : return NULL;
437 : }
438 1003961 : for(uint64_t i = 1; i <= max && i; ++i){
439 1003942 : buf[i] = 1;
440 : }
441 19 : if(!nut_u64_make_factorial_tbl(k, modulus, bits, max_denom, factorials, inv_factorials)){
442 : return NULL;
443 : }
444 1003942 : for(uint64_t n = 2; n <= max && n; ++n){//TODO: optimize loop
445 1003923 : if(nut_Bitarray_get(is_c_buf, n)){
446 924579 : continue;
447 : }
448 2940414 : for(uint64_t m = n; m <= max;){
449 2861070 : nut_Bitarray_set(is_c_buf, m, 1);
450 2861070 : uint64_t m1 = m, a = 0;
451 6497805 : while(m1%n == 0){
452 3636735 : m1 /= n;
453 3636735 : ++a;
454 : }
455 2861070 : uint64_t term = buf[m]*factorials[a + k - 1]%modulus;
456 2861070 : term = term*inv_factorials[k - 1]%modulus;
457 2861070 : buf[m] = term*inv_factorials[a]%modulus;
458 2861070 : if(__builtin_add_overflow(m, n, &m)){
459 : break;
460 : }
461 : }
462 : }
463 : return buf;
464 : }
465 :
466 6 : static uint64_t *sieve_dk_u64(uint64_t max, uint64_t k){
467 6 : uint64_t *buf = malloc((max + 1)*sizeof(uint64_t));
468 12 : uint8_t *is_c_buf [[gnu::cleanup(cleanup_free)]] = calloc(max/8 + 1, sizeof(uint8_t));
469 6 : if(!buf || !is_c_buf){
470 : return NULL;
471 : }
472 1005006 : for(uint64_t i = 1; i <= max && i; ++i){
473 1005000 : buf[i] = 1;
474 : }
475 1005000 : for(uint64_t n = 2; n <= max && n; ++n){//TODO: optimize loop
476 1004994 : if(nut_Bitarray_get(is_c_buf, n)){
477 925656 : continue;
478 : }
479 2943676 : for(uint64_t m = n; m <= max;){
480 2864338 : nut_Bitarray_set(is_c_buf, m, 1);
481 2864338 : uint64_t m1 = m, a = 0;
482 6505342 : while(m1%n == 0){
483 3641004 : m1 /= n;
484 3641004 : ++a;
485 : }
486 2864338 : buf[m] *= nut_u64_binom(a + k - 1, k - 1);
487 2864338 : if(__builtin_add_overflow(m, n, &m)){
488 : break;
489 : }
490 : }
491 : }
492 : return buf;
493 : }
494 :
495 25 : uint64_t *nut_sieve_dk(uint64_t max, uint64_t k, uint64_t modulus){
496 25 : return modulus ? sieve_dk_mod(max, k, modulus) : sieve_dk_u64(max, k);
497 : }
498 :
499 2 : uint64_t *nut_sieve_phi(uint64_t max){
500 2 : uint64_t *buf = malloc((max + 1)*sizeof(uint64_t));
501 2 : if(!buf){
502 : return NULL;
503 : }
504 1001002 : for(uint64_t i = 1; i <= max && i; ++i){
505 1001000 : buf[i] = 1;
506 : }
507 1001000 : for(uint64_t n = 2; n <= max && n; ++n){//TODO: optimize loop
508 1000998 : if(buf[n] != 1){
509 922332 : continue;
510 : }
511 2934500 : for(uint64_t m = n; m <= max;){
512 2855834 : uint64_t a = n*n;
513 3629496 : while(m%a == 0){
514 773662 : a *= n;
515 : }
516 2855834 : buf[m] *= a/n - a/n/n;
517 2855834 : if(__builtin_add_overflow(m, n, &m)){
518 : break;
519 : }
520 : }
521 : }
522 : return buf;
523 : }
524 :
525 2 : uint64_t *nut_sieve_carmichael(uint64_t max){
526 2 : uint64_t *buf = malloc((max + 1)*sizeof(uint64_t));
527 2 : if(!buf){
528 : return NULL;
529 : }
530 1010002 : for(uint64_t i = 1; i <= max && i; ++i){
531 1010000 : buf[i] = 1;
532 : }
533 : uint64_t t = 4, b = 1;
534 34 : while((t >> 1) <= max){
535 505032 : for(uint64_t m = t >> 1; m <= max; m += t){
536 505000 : buf[m] = b;
537 : }
538 32 : t <<= 1;
539 32 : if(t != 16){
540 30 : b <<= 1;
541 : }
542 : }
543 505000 : for(uint64_t n = 1; !__builtin_add_overflow(n, 2, &n) && n <= max;){//TODO: optimize loop
544 504998 : if(buf[n] != 1){
545 425273 : continue;
546 : }
547 2452733 : for(uint64_t k = 1, m = n; k <= max/n; ++k, m += n){
548 2373008 : uint64_t a = n - 1;
549 2648616 : for(uint64_t j = k; j%n == 0; j /= n){
550 275608 : a *= n;
551 : }
552 2373008 : buf[m] = nut_i64_lcm(buf[m], a);
553 : }
554 : }
555 : return buf;
556 : }
557 :
558 3 : uint8_t *nut_sieve_mobius(uint64_t max){
559 3 : uint64_t buf_len = max/4 + 1;
560 3 : uint8_t *buf = malloc(buf_len*sizeof(uint64_t));
561 3 : if(!buf){
562 : return NULL;
563 : }
564 : // Initialize buf to be 0 for 0, 1 for 1, and 2 for every other n <= max.
565 : // 2 is not a valid output value, so we use it here to mark a number as prime
566 : // until they are marked off by being the factor of some other number.
567 3 : buf[0] = 0xA4u;
568 3 : memset(buf + 1, 0xAA, (buf_len - 1)*sizeof(uint8_t));
569 2160428 : for(uint64_t n = 2; n <= max && n; ++n){
570 2160425 : if(((buf[n/4] >> (n%4*2)) & 3) != 2){
571 1991858 : continue;
572 : }
573 168567 : buf[n/4] ^= 1ull << (n%4*2);
574 168567 : uint64_t m;
575 168567 : if(__builtin_mul_overflow(n, 2, &m)){
576 0 : continue;
577 : }
578 : // First we visit all multiples of n and flip their mobius value
579 : // We need to send 2 to 3, 1 to 3, 3 to 1, and 0 to 0.
580 : // We can do this by using an xor mask so we only have to
581 : // get x from the array and then xor the array with this mask.
582 6178061 : while(m <= max){
583 6009494 : uint8_t x = (buf[m/4] >> (m%4*2)) & 3;
584 6009494 : x = (0x98u >> (x*2)) & 3;// this constant acts as a packed lookup table for the xor mask for x
585 6009494 : buf[m/4] ^= x << (m%4*2);
586 6009494 : if(__builtin_add_overflow(m, n, &m)){
587 : break;
588 : }
589 : }
590 168567 : uint64_t n2;
591 : // Here we visit all multiples of n**2 and set their mobius value to 0,
592 : // since they are not square free.
593 168567 : if(!__builtin_mul_overflow(n, n, &n2)){
594 1145197 : for(uint64_t m = n2; m <= max;){
595 976630 : buf[m/4] &= ~(3u << (m%4*2));
596 976630 : if(__builtin_add_overflow(m, n2, &m)){
597 : break;
598 : }
599 : }
600 : }
601 : }
602 : return buf;
603 : }
604 :
605 1 : int64_t *nut_compute_mertens_range(uint64_t max, const uint8_t mobius[static max/4 + 1]){
606 1 : int64_t *buf = malloc((max + 1)*sizeof(int64_t));
607 1 : if(!buf){
608 : return NULL;
609 : }
610 1 : buf[0] = 0;
611 1000001 : for(uint64_t n = 1; n <= max; ++n){
612 1000000 : int64_t v = nut_Bitfield2_arr_get(mobius, n);
613 1000000 : if(v == 3){
614 303857 : v = -1;
615 : }
616 1000000 : buf[n] = buf[n - 1] + v;
617 : }
618 : return buf;
619 : }
620 :
621 518 : static void mark_is_composite(uint8_t *is_composite, uint64_t q, uint64_t r, uint64_t q_ub){
622 : /* We need to find the first multiple of the prime 30*q + r which could possibly be prime and is at least the prime squared.
623 : * Since the prime squared is still coprime to the wheel size, we start at the prime squared.
624 : * (30*q + r)**2 = 30*30*q**2 + 30*2*q*r + r**2 = 30*(30*q**2 + 2*q*r) + r**2
625 : */
626 518 : uint64_t mq;
627 518 : switch(r){
628 55 : case 1:
629 55 : mq = q*(30*q + 2) + 0;// mr = 1
630 8437 : for(; mq + 28*q + 0 < q_ub; mq += 30*q + 1){
631 8382 : is_composite[mq] |= 0x1;// mr = 1
632 8382 : is_composite[mq + 6*q + 0] |= 0x2;// mr = 7
633 8382 : is_composite[mq + 10*q + 0] |= 0x4;// mr = 11
634 8382 : is_composite[mq + 12*q + 0] |= 0x8;// mr = 13
635 8382 : is_composite[mq + 16*q + 0] |= 0x10;// mr = 17
636 8382 : is_composite[mq + 18*q + 0] |= 0x20;// mr = 19
637 8382 : is_composite[mq + 22*q + 0] |= 0x40;// mr = 23
638 8382 : is_composite[mq + 28*q + 0] |= 0x80;// mr = 29
639 : }
640 55 : if(mq >= q_ub){break;}
641 49 : is_composite[mq] |= 0x1;// mr = 1
642 49 : if((mq += 6*q + 0) >= q_ub){break;}
643 39 : is_composite[mq] |= 0x2;// mr = 7
644 39 : if((mq += 4*q + 0) >= q_ub){break;}
645 30 : is_composite[mq] |= 0x4;// mr = 11
646 30 : if((mq += 2*q + 0) >= q_ub){break;}
647 21 : is_composite[mq] |= 0x8;// mr = 13
648 21 : if((mq += 4*q + 0) >= q_ub){break;}
649 15 : is_composite[mq] |= 0x10;// mr = 17
650 15 : if((mq += 2*q + 0) >= q_ub){break;}
651 12 : is_composite[mq] |= 0x20;// mr = 19
652 12 : if((mq += 4*q + 0) >= q_ub){break;}
653 3 : is_composite[mq] |= 0x40;// mr = 23
654 3 : mq += 6*q + 0; break; // If we get here, we have mq >= q_ub by the fact that we made it out of the for loop
655 74 : case 7:
656 74 : mq = q*(30*q + 14) + 1;// mr = 19
657 23187 : for(; mq + 24*q + 6 < q_ub; mq += 30*q + 7){
658 23113 : is_composite[mq] |= 0x20;// mr = 19
659 23113 : is_composite[mq + 4*q + 1] |= 0x10;// mr = 17
660 23113 : is_composite[mq + 6*q + 2] |= 0x1;// mr = 1
661 23113 : is_composite[mq + 10*q + 2] |= 0x80;// mr = 29
662 23113 : is_composite[mq + 12*q + 3] |= 0x8;// mr = 13
663 23113 : is_composite[mq + 16*q + 4] |= 0x4;// mr = 11
664 23113 : is_composite[mq + 22*q + 5] |= 0x40;// mr = 23
665 23113 : is_composite[mq + 24*q + 6] |= 0x2;// mr = 7
666 : }
667 74 : if(mq >= q_ub){break;}
668 70 : is_composite[mq] |= 0x20;// mr = 19
669 70 : if((mq += 4*q + 1) >= q_ub){break;}
670 61 : is_composite[mq] |= 0x10;// mr = 17
671 61 : if((mq += 2*q + 1) >= q_ub){break;}
672 61 : is_composite[mq] |= 0x1;// mr = 1
673 61 : if((mq += 4*q + 0) >= q_ub){break;}
674 46 : is_composite[mq] |= 0x80;// mr = 29
675 46 : if((mq += 2*q + 1) >= q_ub){break;}
676 40 : is_composite[mq] |= 0x8;// mr = 13
677 40 : if((mq += 4*q + 1) >= q_ub){break;}
678 34 : is_composite[mq] |= 0x4;// mr = 11
679 34 : if((mq += 6*q + 1) >= q_ub){break;}
680 6 : is_composite[mq] |= 0x40;// mr = 23
681 6 : mq += 2*q + 1; break; // If we get here, we have mq >= q_ub by the fact that we made it out of the for loop
682 68 : case 11:
683 68 : mq = q*(30*q + 22) + 4;// mr = 1
684 17386 : for(; mq + 26*q + 9 < q_ub; mq += 30*q + 11){
685 17318 : is_composite[mq] |= 0x1;// mr = 1
686 17318 : is_composite[mq + 2*q + 0] |= 0x40;// mr = 23
687 17318 : is_composite[mq + 6*q + 2] |= 0x2;// mr = 7
688 17318 : is_composite[mq + 8*q + 2] |= 0x80;// mr = 29
689 17318 : is_composite[mq + 12*q + 4] |= 0x8;// mr = 13
690 17318 : is_composite[mq + 18*q + 6] |= 0x20;// mr = 19
691 17318 : is_composite[mq + 20*q + 7] |= 0x4;// mr = 11
692 17318 : is_composite[mq + 26*q + 9] |= 0x10;// mr = 17
693 : }
694 68 : if(mq >= q_ub){break;}
695 55 : is_composite[mq] |= 0x1;// mr = 1
696 55 : if((mq += 2*q + 0) >= q_ub){break;}
697 52 : is_composite[mq] |= 0x40;// mr = 23
698 52 : if((mq += 4*q + 2) >= q_ub){break;}
699 40 : is_composite[mq] |= 0x2;// mr = 7
700 40 : if((mq += 2*q + 0) >= q_ub){break;}
701 31 : is_composite[mq] |= 0x80;// mr = 29
702 31 : if((mq += 4*q + 2) >= q_ub){break;}
703 31 : is_composite[mq] |= 0x8;// mr = 13
704 31 : if((mq += 6*q + 2) >= q_ub){break;}
705 19 : is_composite[mq] |= 0x20;// mr = 19
706 19 : if((mq += 2*q + 1) >= q_ub){break;}
707 10 : is_composite[mq] |= 0x4;// mr = 11
708 10 : mq += 6*q + 2; break; // If we get here, we have mq >= q_ub by the fact that we made it out of the for loop
709 62 : case 13:
710 62 : mq = q*(30*q + 26) + 5;// mr = 19
711 15679 : for(; mq + 28*q + 12 < q_ub; mq += 30*q + 13){
712 15617 : is_composite[mq] |= 0x20;// mr = 19
713 15617 : is_composite[mq + 4*q + 2] |= 0x4;// mr = 11
714 15617 : is_composite[mq + 6*q + 3] |= 0x2;// mr = 7
715 15617 : is_composite[mq + 10*q + 4] |= 0x80;// mr = 29
716 15617 : is_composite[mq + 16*q + 7] |= 0x10;// mr = 17
717 15617 : is_composite[mq + 18*q + 8] |= 0x8;// mr = 13
718 15617 : is_composite[mq + 24*q + 11] |= 0x1;// mr = 1
719 15617 : is_composite[mq + 28*q + 12] |= 0x40;// mr = 23
720 : }
721 62 : if(mq >= q_ub){break;}
722 55 : is_composite[mq] |= 0x20;// mr = 19
723 55 : if((mq += 4*q + 2) >= q_ub){break;}
724 37 : is_composite[mq] |= 0x4;// mr = 11
725 37 : if((mq += 2*q + 1) >= q_ub){break;}
726 30 : is_composite[mq] |= 0x2;// mr = 7
727 30 : if((mq += 4*q + 1) >= q_ub){break;}
728 21 : is_composite[mq] |= 0x80;// mr = 29
729 21 : if((mq += 6*q + 3) >= q_ub){break;}
730 15 : is_composite[mq] |= 0x10;// mr = 17
731 15 : if((mq += 2*q + 1) >= q_ub){break;}
732 9 : is_composite[mq] |= 0x8;// mr = 13
733 9 : if((mq += 6*q + 3) >= q_ub){break;}
734 0 : is_composite[mq] |= 0x1;// mr = 1
735 0 : mq += 4*q + 1; break; // If we get here, we have mq >= q_ub by the fact that we made it out of the for loop
736 68 : case 17:
737 68 : mq = q*(30*q + 34) + 9;// mr = 19
738 12903 : for(; mq + 26*q + 15 < q_ub; mq += 30*q + 17){
739 12835 : is_composite[mq] |= 0x20;// mr = 19
740 12835 : is_composite[mq + 2*q + 1] |= 0x40;// mr = 23
741 12835 : is_composite[mq + 6*q + 4] |= 0x1;// mr = 1
742 12835 : is_composite[mq + 12*q + 7] |= 0x8;// mr = 13
743 12835 : is_composite[mq + 14*q + 8] |= 0x10;// mr = 17
744 12835 : is_composite[mq + 20*q + 11] |= 0x80;// mr = 29
745 12835 : is_composite[mq + 24*q + 14] |= 0x2;// mr = 7
746 12835 : is_composite[mq + 26*q + 15] |= 0x4;// mr = 11
747 : }
748 68 : if(mq >= q_ub){break;}
749 61 : is_composite[mq] |= 0x20;// mr = 19
750 61 : if((mq += 2*q + 1) >= q_ub){break;}
751 55 : is_composite[mq] |= 0x40;// mr = 23
752 55 : if((mq += 4*q + 3) >= q_ub){break;}
753 52 : is_composite[mq] |= 0x1;// mr = 1
754 52 : if((mq += 6*q + 3) >= q_ub){break;}
755 34 : is_composite[mq] |= 0x8;// mr = 13
756 34 : if((mq += 2*q + 1) >= q_ub){break;}
757 30 : is_composite[mq] |= 0x10;// mr = 17
758 30 : if((mq += 6*q + 3) >= q_ub){break;}
759 15 : is_composite[mq] |= 0x80;// mr = 29
760 15 : if((mq += 4*q + 3) >= q_ub){break;}
761 3 : is_composite[mq] |= 0x2;// mr = 7
762 3 : mq += 2*q + 1; break; // If we get here, we have mq >= q_ub by the fact that we made it out of the for loop
763 58 : case 19:
764 58 : mq = q*(30*q + 38) + 12;// mr = 1
765 10445 : for(; mq + 28*q + 17 < q_ub; mq += 30*q + 19){
766 10387 : is_composite[mq] |= 0x1;// mr = 1
767 10387 : is_composite[mq + 4*q + 2] |= 0x10;// mr = 17
768 10387 : is_composite[mq + 10*q + 6] |= 0x4;// mr = 11
769 10387 : is_composite[mq + 12*q + 7] |= 0x20;// mr = 19
770 10387 : is_composite[mq + 18*q + 11] |= 0x8;// mr = 13
771 10387 : is_composite[mq + 22*q + 13] |= 0x80;// mr = 29
772 10387 : is_composite[mq + 24*q + 15] |= 0x2;// mr = 7
773 10387 : is_composite[mq + 28*q + 17] |= 0x40;// mr = 23
774 : }
775 58 : if(mq >= q_ub){break;}
776 55 : is_composite[mq] |= 0x1;// mr = 1
777 55 : if((mq += 4*q + 2) >= q_ub){break;}
778 55 : is_composite[mq] |= 0x10;// mr = 17
779 55 : if((mq += 6*q + 4) >= q_ub){break;}
780 30 : is_composite[mq] |= 0x4;// mr = 11
781 30 : if((mq += 2*q + 1) >= q_ub){break;}
782 27 : is_composite[mq] |= 0x20;// mr = 19
783 27 : if((mq += 6*q + 4) >= q_ub){break;}
784 21 : is_composite[mq] |= 0x8;// mr = 13
785 21 : if((mq += 4*q + 2) >= q_ub){break;}
786 15 : is_composite[mq] |= 0x80;// mr = 29
787 15 : if((mq += 2*q + 2) >= q_ub){break;}
788 12 : is_composite[mq] |= 0x2;// mr = 7
789 12 : mq += 4*q + 2; break; // If we get here, we have mq >= q_ub by the fact that we made it out of the for loop
790 68 : case 23:
791 68 : mq = q*(30*q + 46) + 17;// mr = 19
792 11330 : for(; mq + 26*q + 20 < q_ub; mq += 30*q + 23){
793 11262 : is_composite[mq] |= 0x20;// mr = 19
794 11262 : is_composite[mq + 6*q + 5] |= 0x2;// mr = 7
795 11262 : is_composite[mq + 8*q + 6] |= 0x40;// mr = 23
796 11262 : is_composite[mq + 14*q + 11] |= 0x4;// mr = 11
797 11262 : is_composite[mq + 18*q + 14] |= 0x8;// mr = 13
798 11262 : is_composite[mq + 20*q + 15] |= 0x80;// mr = 29
799 11262 : is_composite[mq + 24*q + 19] |= 0x1;// mr = 1
800 11262 : is_composite[mq + 26*q + 20] |= 0x10;// mr = 17
801 : }
802 68 : if(mq >= q_ub){break;}
803 58 : is_composite[mq] |= 0x20;// mr = 19
804 58 : if((mq += 6*q + 5) >= q_ub){break;}
805 43 : is_composite[mq] |= 0x2;// mr = 7
806 43 : if((mq += 2*q + 1) >= q_ub){break;}
807 37 : is_composite[mq] |= 0x40;// mr = 23
808 37 : if((mq += 6*q + 5) >= q_ub){break;}
809 19 : is_composite[mq] |= 0x4;// mr = 11
810 19 : if((mq += 4*q + 3) >= q_ub){break;}
811 10 : is_composite[mq] |= 0x8;// mr = 13
812 10 : if((mq += 2*q + 1) >= q_ub){break;}
813 7 : is_composite[mq] |= 0x80;// mr = 29
814 7 : if((mq += 4*q + 4) >= q_ub){break;}
815 6 : is_composite[mq] |= 0x1;// mr = 1
816 6 : mq += 2*q + 1; break; // If we get here, we have mq >= q_ub by the fact that we made it out of the for loop
817 65 : case 29:
818 65 : mq = q*(30*q + 58) + 28;// mr = 1
819 9839 : for(; mq + 24*q + 23 < q_ub; mq += 30*q + 29){
820 9774 : is_composite[mq] |= 0x1;// mr = 1
821 9774 : is_composite[mq + 2*q + 1] |= 0x80;// mr = 29
822 9774 : is_composite[mq + 8*q + 7] |= 0x40;// mr = 23
823 9774 : is_composite[mq + 12*q + 11] |= 0x20;// mr = 19
824 9774 : is_composite[mq + 14*q + 13] |= 0x10;// mr = 17
825 9774 : is_composite[mq + 18*q + 17] |= 0x8;// mr = 13
826 9774 : is_composite[mq + 20*q + 19] |= 0x4;// mr = 11
827 9774 : is_composite[mq + 24*q + 23] |= 0x2;// mr = 7
828 : }
829 65 : if(mq >= q_ub){break;}
830 52 : is_composite[mq] |= 0x1;// mr = 1
831 52 : if((mq += 2*q + 1) >= q_ub){break;}
832 49 : is_composite[mq] |= 0x80;// mr = 29
833 49 : if((mq += 6*q + 6) >= q_ub){break;}
834 42 : is_composite[mq] |= 0x40;// mr = 23
835 42 : if((mq += 4*q + 4) >= q_ub){break;}
836 36 : is_composite[mq] |= 0x20;// mr = 19
837 36 : if((mq += 2*q + 2) >= q_ub){break;}
838 36 : is_composite[mq] |= 0x10;// mr = 17
839 36 : if((mq += 4*q + 4) >= q_ub){break;}
840 21 : is_composite[mq] |= 0x8;// mr = 13
841 21 : if((mq += 2*q + 2) >= q_ub){break;}
842 12 : is_composite[mq] |= 0x4;// mr = 11
843 12 : mq += 4*q + 4; break; // If we get here, we have mq >= q_ub by the fact that we made it out of the for loop
844 0 : default:
845 0 : __builtin_unreachable();
846 : }
847 518 : }
848 :
849 4 : uint8_t *nut_sieve_is_composite(uint64_t max){
850 4 : uint64_t is_composite_len = max/30 + 1;
851 4 : uint8_t *is_composite = calloc(is_composite_len, sizeof(uint64_t));
852 4 : if(!is_composite){
853 : return NULL;
854 : }
855 4 : uint64_t q_ub = is_composite_len;
856 4 : uint64_t p_max = nut_u64_nth_root(max, 2);
857 4 : uint64_t q_max = (p_max - 1)/30;
858 4 : uint64_t q = 0, r = 7;
859 108 : while(q <= q_max){
860 104 : uint8_t flags = is_composite[q];
861 104 : switch(r){
862 100 : case 1:
863 100 : if(!(flags & 0x1)){mark_is_composite(is_composite, q, r, q_ub);}
864 104 : r = 7; [[fallthrough]];
865 104 : case 7:
866 104 : if(!(flags & 0x2)){mark_is_composite(is_composite, q, r, q_ub);}
867 104 : r = 11; [[fallthrough]];
868 104 : case 11:
869 104 : if(!(flags & 0x4)){mark_is_composite(is_composite, q, r, q_ub);}
870 104 : r = 13; [[fallthrough]];
871 104 : case 13:
872 104 : if(!(flags & 0x8)){mark_is_composite(is_composite, q, r, q_ub);}
873 104 : r = 17; [[fallthrough]];
874 104 : case 17:
875 104 : if(!(flags & 0x10)){mark_is_composite(is_composite, q, r, q_ub);}
876 104 : r = 19; [[fallthrough]];
877 104 : case 19:
878 104 : if(!(flags & 0x20)){mark_is_composite(is_composite, q, r, q_ub);}
879 104 : r = 23; [[fallthrough]];
880 104 : case 23:
881 104 : if(!(flags & 0x40)){mark_is_composite(is_composite, q, r, q_ub);}
882 104 : r = 29; [[fallthrough]];
883 104 : case 29:
884 104 : if(!(flags & 0x80)){mark_is_composite(is_composite, q, r, q_ub);}
885 : r = 1;
886 : ++q;
887 : break;
888 0 : default:
889 0 : __builtin_unreachable();
890 : }
891 : }
892 : return is_composite;
893 : }
894 :
895 1000000 : bool nut_is_composite(uint64_t n, const uint8_t buf[static n/30 + 1]){
896 1000000 : if(n < 6){
897 5 : return n != 2 && n != 3 && n != 5;
898 : }
899 999995 : uint64_t r = n%30;
900 999995 : uint64_t q = n/30;
901 999995 : switch(r){
902 33333 : case 1: return buf[q] & 0x1;
903 33334 : case 7: return buf[q] & 0x2;
904 33333 : case 11: return buf[q] & 0x4;
905 33333 : case 13: return buf[q] & 0x8;
906 33333 : case 17: return buf[q] & 0x10;
907 33333 : case 19: return buf[q] & 0x20;
908 33333 : case 23: return buf[q] & 0x40;
909 33333 : case 29: return buf[q] & 0x80;
910 : default: return true;
911 : }
912 : }
913 :
914 2 : uint64_t *nut_compute_pi_range(uint64_t max, const uint8_t buf[static max/30 + 1]){
915 : /* buf[i] tells us the primality of 30*i + 1, 30*i + 7, 30*i + 11, ..., 30*i + 29.
916 : * Thus, we want res[i] to tell us the number of prime numbers from 1 up to 30*(i + 1)
917 : * and so in general res[i] = popcount(buf[i]) + res[i - 1].
918 : * However, res[0], the number of primes up to 30, is instead 10.
919 : * Finally, since max is inclusive, we must have res[max/30 - 1].
920 : * The -1 is because since res[0] is the number of primes up to 30, for arguments i < 30 we have to
921 : * special case their lookup, so the domain essentially starts at 30.
922 : */
923 2 : uint64_t res_len = max/30;
924 2 : uint64_t *res = calloc(res_len, sizeof(uint64_t));
925 2 : if(!res){
926 : return NULL;
927 : }
928 2 : res[0] = 10;
929 66666 : for(uint64_t i = 1, acc = 10; i < res_len; ++i){
930 66664 : res[i] = acc += __builtin_popcount(0xFF&~buf[i]);
931 : }
932 : return res;
933 : }
934 :
935 1002000 : uint64_t nut_compute_pi_from_tables(uint64_t n, const uint64_t pi_table[restrict static n/30], const uint8_t buf[restrict static n/30 + 1]){
936 1002000 : if(n < 30){
937 : uint64_t res;
938 400 : for(res = 0; nut_small_primes[res] <= n; ++res);
939 : return res;
940 : }
941 1001942 : uint64_t q = n/30;
942 1001942 : uint64_t r = n%30;
943 1001942 : uint8_t mask;
944 1001942 : if(!r){mask = 0x00;}
945 968545 : else if(r < 7){mask = 0x01;}
946 768148 : else if(r < 11){mask = 0x03;}
947 634545 : else if(r < 13){mask = 0x07;}
948 567752 : else if(r < 17){mask = 0x0f;}
949 434161 : else if(r < 19){mask = 0x1f;}
950 367355 : else if(r < 23){mask = 0x3f;}
951 233767 : else if(r < 29){mask = 0x7f;}
952 33391 : else{mask = 0xff;}
953 1001942 : return pi_table[q-1] + __builtin_popcount(mask&~buf[q]);
954 : }
955 :
956 21 : static uint64_t *copy_small_primes(uint64_t max, uint64_t *_num_primes){
957 21 : uint64_t n = 0;
958 139 : while(n < 25 && nut_small_primes[n] <= max){
959 118 : ++n;
960 : }
961 21 : uint64_t *primes = malloc(n*sizeof(uint64_t));
962 21 : if(!primes){
963 : return NULL;
964 : }
965 21 : memcpy(primes, nut_small_primes, n*sizeof(uint64_t));
966 21 : *_num_primes = n;
967 21 : return primes;
968 : }
969 :
970 23 : uint64_t *nut_sieve_primes(uint64_t max, uint64_t *_num_primes){
971 23 : if(max <= 100){
972 21 : return copy_small_primes(max, _num_primes);
973 : }
974 2 : uint64_t *primes = malloc((size_t)nut_max_primes_le(max)*sizeof(uint64_t));
975 2 : if(!primes){
976 : return NULL;
977 : }
978 2 : uint8_t *is_composite = nut_sieve_is_composite(max);
979 2 : if(!is_composite){
980 0 : free(primes);
981 0 : return NULL;
982 : }
983 2 : memcpy(primes, nut_small_primes, 3*sizeof(uint64_t));
984 2 : uint64_t num_primes = 3;
985 2 : uint64_t q = 0;
986 2 : uint64_t r = 7;
987 33368 : while(30*q + 29 <= max){
988 33366 : uint8_t flags = is_composite[q];
989 33366 : switch(r){
990 33364 : case 1:
991 33364 : if(!(flags & 0x1)){primes[num_primes++] = 30*q + r;}
992 33366 : r = 7; [[fallthrough]];
993 33366 : case 7:
994 33366 : if(!(flags & 0x2)){primes[num_primes++] = 30*q + r;}
995 33366 : r = 11; [[fallthrough]];
996 33366 : case 11:
997 33366 : if(!(flags & 0x4)){primes[num_primes++] = 30*q + r;}
998 33366 : r = 13; [[fallthrough]];
999 33366 : case 13:
1000 33366 : if(!(flags & 0x8)){primes[num_primes++] = 30*q + r;}
1001 33366 : r = 17; [[fallthrough]];
1002 33366 : case 17:
1003 33366 : if(!(flags & 0x10)){primes[num_primes++] = 30*q + r;}
1004 33366 : r = 19; [[fallthrough]];
1005 33366 : case 19:
1006 33366 : if(!(flags & 0x20)){primes[num_primes++] = 30*q + r;}
1007 33366 : r = 23; [[fallthrough]];
1008 33366 : case 23:
1009 33366 : if(!(flags & 0x40)){primes[num_primes++] = 30*q + r;}
1010 33366 : r = 29; [[fallthrough]];
1011 33366 : case 29:
1012 33366 : if(!(flags & 0x80)){primes[num_primes++] = 30*q + r;}
1013 : r = 1;
1014 : ++q;
1015 : break;
1016 0 : default:
1017 0 : __builtin_unreachable();
1018 : }
1019 : }
1020 4 : do{
1021 2 : if(30*q + 1 > max){break;}
1022 2 : uint8_t flags = is_composite[q];
1023 2 : if(!(flags & 0x1)){primes[num_primes++] = 30*q + 1;}
1024 2 : if(30*q + 7 > max){break;}
1025 2 : if(!(flags & 0x2)){primes[num_primes++] = 30*q + 7;}
1026 2 : if(30*q + 11 > max){break;}
1027 0 : if(!(flags & 0x4)){primes[num_primes++] = 30*q + 11;}
1028 0 : if(30*q + 13 > max){break;}
1029 0 : if(!(flags & 0x8)){primes[num_primes++] = 30*q + 13;}
1030 0 : if(30*q + 17 > max){break;}
1031 0 : if(!(flags & 0x10)){primes[num_primes++] = 30*q + 17;}
1032 0 : if(30*q + 19 > max){break;}
1033 0 : if(!(flags & 0x20)){primes[num_primes++] = 30*q + 19;}
1034 0 : if(30*q + 23 > max){break;}
1035 0 : if(!(flags & 0x40)){primes[num_primes++] = 30*q + 23;}
1036 : if(30*q + 29 > max){break;}
1037 : if(!(flags & 0x80)){primes[num_primes++] = 30*q + 29;}
1038 : }while(0);
1039 2 : *_num_primes = num_primes;
1040 2 : free(is_composite);
1041 2 : return realloc(primes, num_primes*sizeof(uint64_t)) ?: primes;
1042 : }
1043 :
|
