File Coverage

sieve.h
Criterion Covered Total %
statement 12 13 92.3
branch 5 14 35.7
condition n/a
subroutine n/a
pod n/a
total 17 27 62.9


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