Line data Source code
1 : #include <inttypes.h>
2 : #include <stdlib.h>
3 : #include <stddef.h>
4 :
5 : #include <nut/modular_math.h>
6 : #include <nut/factorization.h>
7 : #include <nut/sieves.h>
8 : #include <nut/segsieves.h>
9 :
10 1 : bool nut_Segsieve_init(nut_Segsieve *self, uint64_t max, uint64_t preferred_bucket_size){
11 1 : self->max = max;
12 1 : self->sqrt_max = nut_u64_nth_root(max, 2);
13 1 : self->preferred_bucket_size = preferred_bucket_size ?: self->sqrt_max;
14 1 : return (self->primes = nut_sieve_primes(self->sqrt_max, &self->num_primes));
15 : }
16 :
17 1 : void nut_Segsieve_destroy(nut_Segsieve *self){
18 1 : free(self->primes);
19 1 : }
20 :
21 1001 : void nut_Segsieve_factorizations(const nut_Segsieve *restrict self, uint64_t a, uint64_t b, size_t pitch, void *buffer){
22 1001000 : for(uint64_t n = a; n < b; ++n){
23 999999 : nut_Factors *fxn = nut_Pitcharr_get(buffer, pitch, n - a);
24 999999 : fxn->num_primes = 1;
25 999999 : fxn->factors[0].prime = n;
26 999999 : fxn->factors[0].power = 1;
27 : }
28 169169 : for(uint64_t i = 0; i < self->num_primes; ++i){
29 168168 : uint64_t p = self->primes[i];
30 2366175 : for(uint64_t m = (a + p - 1)/p*p; m < b; m += p){
31 2198007 : nut_Factors *fxn = nut_Pitcharr_get(buffer, pitch, m - a);
32 2198007 : fxn->factors[fxn->num_primes].prime = p;
33 2198007 : fxn->factors[fxn->num_primes].power = 1;
34 2198007 : fxn->factors[0].prime /= p;
35 2970918 : while(fxn->factors[0].prime%p == 0){
36 772911 : fxn->factors[0].prime /= p;
37 772911 : fxn->factors[fxn->num_primes].power++;
38 : }
39 2198007 : fxn->num_primes++;
40 : }
41 : }
42 1001 : }
43 :
44 2 : void *nut_Segsieve_factorizations_mkbuffer(const nut_Segsieve *self, size_t *pitch){
45 2 : if(!*pitch){
46 2 : *pitch = offsetof(nut_Factors, factors) + (nut_max_prime_divs(self->max) + 1)*sizeof(uint64_t)*2;
47 : }
48 2 : return malloc(*pitch*self->preferred_bucket_size);
49 : }
50 :
|
