1 #ifndef MPU_SIEVE_H
2 #define MPU_SIEVE_H
3 
4 #include "ptypes.h"
5 #define FUNC_ctz 1
6 #include "util.h"
7 
8 extern unsigned char* sieve_erat30(UV end);
9 extern int sieve_segment_partial(unsigned char* mem, UV startd, UV endd, UV depth);
10 extern int sieve_segment(unsigned char* mem, UV startd, UV endd);
11 extern void* start_segment_primes(UV low, UV high, unsigned char** segmentmem);
12 extern int next_segment_primes(void* vctx, UV* base, UV* low, UV* high);
13 extern void end_segment_primes(void* vctx);
14 
15 extern void* array_of_primes_in_range(UV* count, UV beg, UV end);
16 
17 static const UV wheel30[] = {1, 7, 11, 13, 17, 19, 23, 29};
18 /* Used for moving between primes */
19 static const unsigned char nextwheel30[30] = {
20     1,  7,  7,  7,  7,  7,  7, 11, 11, 11, 11, 13, 13, 17, 17,
21    17, 17, 19, 19, 23, 23, 23, 23, 29, 29, 29, 29, 29, 29,  1 };
22 static const unsigned char prevwheel30[30] = {
23    29, 29,  1,  1,  1,  1,  1,  1,  7,  7,  7,  7, 11, 11, 13,
24    13, 13, 13, 17, 17, 19, 19, 19, 19, 23, 23, 23, 23, 23, 23 };
25 /* The bit mask within a byte */
26 static const unsigned char masktab30[30] = {
27     0,  1,  0,  0,  0,  0,  0,  2,  0,  0,  0,  4,  0,  8,  0,
28     0,  0, 16,  0, 32,  0,  0,  0, 64,  0,  0,  0,  0,  0,128  };
29 /* Inverse of masktab30 */
30 static const unsigned char imask30[129] = {
31 0,1,7,0,11,0,0,0,13,0,0,0,0,0,0,0,17,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,19,
32 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,23,
33 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
34 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,29};
35 /* Add this to a number and you'll ensure you're on a wheel location */
36 static const unsigned char distancewheel30[30] =
37     {1,0,5,4,3,2,1,0,3,2,1,0,1,0,3,2,1,0,1,0,3,2,1,0,5,4,3,2,1,0};
38 /* add this to n to get to the next wheel location */
39 static const unsigned char wheeladvance30[30] =
40     {1,6,5,4,3,2,1,4,3,2,1,2,1,4,3,2,1,2,1,4,3,2,1,6,5,4,3,2,1,2};
41 /* subtract this from n to get to the previous wheel location */
42 static const unsigned char wheelretreat30[30] =
43     {1,2,1,2,3,4,5,6,1,2,3,4,1,2,1,2,3,4,1,2,1,2,3,4,1,2,3,4,5,6};
44 /* Given a sieve byte, this indicates the first zero */
45 static const unsigned char nextzero30[256] =
46   {0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,5,0,1,0,2,0,1,
47    0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,6,0,1,0,2,0,1,0,3,0,1,0,2,
48    0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,5,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,
49    0,2,0,1,0,3,0,1,0,2,0,1,0,7,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,
50    0,1,0,2,0,1,0,5,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,
51    0,6,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,5,0,1,0,2,
52    0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,8};
53 /* At this m (p-30*(p/30)), OR with this to clear previous entries */
54 static const unsigned char clearprev30[30] =
55     {  0,  0,  1,  1,  1,  1,  1,  1,  3,  3,  3,  3,  7,  7, 15,
56       15, 15, 15, 31, 31, 63, 63, 63, 63,127,127,127,127,127,127};
57 
58 
59 #ifdef FUNC_is_prime_in_sieve
is_prime_in_sieve(const unsigned char * sieve,UV p)60 static int is_prime_in_sieve(const unsigned char* sieve, UV p) {
61   UV d = p/30;
62   UV m = p - d*30;
63   /* If m isn't part of the wheel, we return 0 */
64   return ( (masktab30[m] != 0) && ((sieve[d] & masktab30[m]) == 0) );
65 }
66 #endif
67 
68 #ifdef FUNC_next_prime_in_sieve
69 /* Will return 0 if it goes past lastp */
next_prime_in_sieve(const unsigned char * sieve,UV p,UV lastp)70 static UV next_prime_in_sieve(const unsigned char* sieve, UV p, UV lastp) {
71   UV d, m;
72   unsigned char s;
73   if (p < 7)
74     return (p < 2) ? 2 : (p < 3) ? 3 : (p < 5) ? 5 : 7;
75   p++;
76   if (p >= lastp) return 0;
77   d = p/30;
78   m = p - d*30;
79   s = sieve[d] | clearprev30[m];
80   while (s == 0xFF) {
81     d++;
82     if (d*30 >= lastp) return 0;
83     s = sieve[d];
84   }
85   return d*30 + wheel30[nextzero30[s]];
86 }
87 #endif
88 #ifdef FUNC_prev_prime_in_sieve
prev_prime_in_sieve(const unsigned char * sieve,UV p)89 static UV prev_prime_in_sieve(const unsigned char* sieve, UV p) {
90   UV d, m;
91   if (p <= 7)
92     return (p <= 2) ? 0 : (p <= 3) ? 2 : (p <= 5) ? 3 : 5;
93   d = p/30;
94   m = p - d*30;
95   do {
96     m = prevwheel30[m];
97     if (m==29) { if (d == 0) return 0;  d--; }
98   } while (sieve[d] & masktab30[m]);
99   return(d*30+m);
100 }
101 #endif
102 
103 #if 0
104 /* Useful macros for the wheel-30 sieve array */
105 #define START_DO_FOR_EACH_SIEVE_PRIME(sieve, base, a, b) \
106   { \
107     const unsigned char* sieve_ = sieve; \
108     UV base_ = base; \
109     UV p = a-base_; \
110     UV l_ = b; \
111     UV d_ = p/30; \
112     UV lastd_ = (l_-base_)/30; \
113     unsigned char bit_, s_ = sieve_[d_] | clearprev30[p-d_*30]; \
114     base_ += d_*30; \
115     while (1) { \
116       if (s_ == 0xFF) { \
117         do { \
118           base_ += 30;  d_++; \
119           if (d_ > lastd_) break; \
120           s_ = sieve_[d_]; \
121         } while (s_ == 0xFF); \
122         if (d_ > lastd_) break; \
123       } \
124       bit_ = nextzero30[s_]; \
125       s_ |= 1 << bit_; \
126       p = base_ + wheel30[bit_]; \
127       if (p > l_ || p < base_) break; /* handle overflow */ \
128       {
129 
130 #define END_DO_FOR_EACH_SIEVE_PRIME \
131       } \
132     } \
133   }
134 #else
135 /* Extract word at a time, good suggestion from Kim Walisch */
136 static const unsigned char wheel240[] = {1,7,11,13,17,19,23,29,31,37,41,43,47,49,53,59,61,67,71,73,77,79,83,89,91,97,101,103,107,109,113,119,121,127,131,133,137,139,143,149,151,157,161,163,167,169,173,179,181,187,191,193,197,199,203,209,211,217,221,223,227,229,233,239};
137 #define START_DO_FOR_EACH_SIEVE_PRIME(sieve, base, a, b) \
138   { \
139     const UV* sieve_ = (const UV*)sieve;        /* word ptr to sieve */ \
140     const UV nperw_ = 30*sizeof(UV); /* nums per word */     \
141     UV base_ = base;                 /* start of sieve n */  \
142     UV b_ = a;                       /* begin value n */     \
143     UV f_ = b;                       /* final value n */     \
144     UV begw_ = (b_-base_)/nperw_;    /* first word */        \
145     UV endw_ = (f_-base_)/nperw_;    /* first word */        \
146     UV sw_, tz_, p; \
147     base_ += begw_*nperw_; \
148     while (begw_ <= endw_) { \
149       sw_ = ~ LEUV(sieve_[begw_]); \
150       while (sw_ != 0) { \
151         tz_ = ctz(sw_); \
152         sw_ &= ~(UVCONST(1) << tz_); \
153         p = base_ + wheel240[tz_]; \
154         if (p > f_) break; \
155         if (p >= b_) {
156 
157 #define END_DO_FOR_EACH_SIEVE_PRIME \
158         } \
159       } \
160       begw_++; \
161       base_ += nperw_; \
162     } \
163   }
164 #endif
165 
166 #define START_DO_FOR_EACH_PRIME(a, b) \
167   { \
168     const unsigned char* sieve_; \
169     UV p  = a; \
170     UV l_ = b; \
171     UV d_ = p/30; \
172     UV lastd_ = l_/30; \
173     unsigned char s_, bit_; \
174     get_prime_cache(l_, &sieve_); \
175     if (p == 2) p = 1; \
176     s_ = sieve_[d_] | clearprev30[p-d_*30]; \
177     while (1) { \
178       if (p < 5) { \
179         p = (p < 2) ? 2 : (p < 3) ? 3 : 5; \
180       } else { \
181         if (s_ == 0xFF) { \
182           do { \
183             d_++; \
184             if (d_ > lastd_) break; \
185             s_ = sieve_[d_]; \
186           } while (s_ == 0xFF); \
187           if (d_ > lastd_) break; \
188         } \
189         bit_ = nextzero30[s_]; \
190         s_ |= 1 << bit_; \
191         p = d_*30 + wheel30[bit_]; \
192         if (p < d_*30) break; \
193       } \
194       if (p > l_) break; \
195       { \
196 
197 #define RETURN_FROM_EACH_PRIME(retstmt) \
198         do { release_prime_cache(sieve_); retstmt; } while (0)
199 
200 #define END_DO_FOR_EACH_PRIME \
201       } \
202     } \
203     release_prime_cache(sieve_); \
204   }
205 
206 #endif
207