chiark / gitweb /
Merge branch 'mdw/rsvr'
[catacomb] / math / primeiter.c
CommitLineData
1d6d3b01
MW
1/* -*-c-*-
2 *
3 * Iterate over small primes efficiently
4 *
5 * (c) 2007 Straylight/Edgeware
6 */
7
8/*----- Licensing notice --------------------------------------------------*
9 *
10 * This file is part of Catacomb.
11 *
12 * Catacomb is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU Library General Public License as
14 * published by the Free Software Foundation; either version 2 of the
15 * License, or (at your option) any later version.
16 *
17 * Catacomb is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Library General Public License for more details.
21 *
22 * You should have received a copy of the GNU Library General Public
23 * License along with Catacomb; if not, write to the Free
24 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
25 * MA 02111-1307, USA.
26 */
27
28/*----- Header files ------------------------------------------------------*/
29
30#include "fibrand.h"
31#include "mp.h"
32#include "pgen.h"
33#include "primeiter.h"
34#include "primetab.h"
35#include "wheel.h"
36
37/*----- Theory ------------------------------------------------------------*
38 *
39 * For small primes, we can just pluck them out of the small primes table.
40 * For larger primes, we can test them individually, or build a sieve or
41 * something, but since we don't know when to stop, that could be tricky.
42 *
43 * We've built a `wheel', as follows. Let %$m$% be the product of the first
44 * %$n$% primes. There are %$\phi(m)$% integers %$n_i$%, with %$0 < n_i <
45 * m$% coprime to %$m$%, and any integer %$j > n$% must be congruent to some
46 * %$n_i$% modulo %$m$%. The wheel itself doesn't list the %$n_i$%, but
47 * rather the differences %$\delta_i = n_i - n_{i-1}$% (wrapping around
48 * appropriately at the ends), so you can just add simple offsets to step
49 * onwards. The wheel assumes you start at 1 and move on round.
50 */
51
52/*----- Main code ---------------------------------------------------------*/
53
54/* --- @wheelsync@ --- *
55 *
56 * Arguments: @primeiter *pi@ = iterator to synchronize
57 * @mp *where@ = value to synchronize
58 *
59 * Returns: ---
60 *
61 * Use: Sets up the wheel index to match the given integer. After
62 * this, we can step along the wheel to find candidate primes.
63 */
64
65static void wheelsync(primeiter *pi, mp *where)
66{
67 mpw w;
68 mp t;
69 mpw rr;
70 mp *r = MP_NEW;
71 unsigned i, n;
72
73 w = WHEELMOD;
74 mp_build(&t, &w, &w + 1);
75 mp_div(0, &r, where, &t);
76 rr = MP_ZEROP(r) ? 0 : r->v[0];
77
78 for (i = 0, n = 1; rr > n; n += wheel[i], i++);
79 w = n - rr;
80 pi->p = mp_add(MP_NEW, where, &t);
81 pi->i = i;
82 pi->r = fibrand_create(0);
83 MP_DROP(r);
84}
85
86/* --- @primeiter_create@ --- *
87 *
88 * Arguments: @primeiter *pi@ = pointer to an iterator structure
89 * @mp *start@ = where to start
90 *
91 * Returns: ---
92 *
93 * Use: Initializes a prime iterator. The first output will be the
94 * smallest prime not less than @start@.
95 */
96
97void primeiter_create(primeiter *pi, mp *start)
98{
99 mpw n;
100 unsigned l, h, m;
101
102 if (!start || MP_CMP(start, <=, MP_TWO))
103 start = MP_TWO;
104
105 if (MP_LEN(start) <= 1) {
106 n = start->v[0];
107 if (n <= MAXPRIME) {
108 l = 0;
109 h = NPRIME;
110 for (;;) {
111 m = l + (h - l)/2;
112 if (primetab[m] == n) break;
113 else if (m == l) { m++; break; }
114 else if (primetab[m] < n) l = m;
115 else h = m;
116 }
117 pi->i = m;
118 pi->mode = PIMODE_PTAB;
119 mp_build(&pi->pp, &pi->w, &pi->w + 1);
120 pi->p = &pi->pp;
121 return;
122 }
123 }
124
125 wheelsync(pi, start);
126 pi->mode = PIMODE_STALL;
127}
128
129/* --- @primeiter_destroy@ --- *
130 *
131 * Arguments: @primeiter *pi@ = pointer to iterator structure
132 *
133 * Returns: ---
134 *
135 * Use: Frees up an iterator structure when it's no longer wanted.
136 */
137
138void primeiter_destroy(primeiter *pi)
139{
140 switch (pi->mode) {
141 case PIMODE_PTAB:
142 break;
143 case PIMODE_STALL:
144 case PIMODE_WHEEL:
145 MP_DROP(pi->p);
146 GR_DESTROY(pi->r);
147 break;
148 default:
149 abort();
150 }
151}
152
153/* --- @primeiter_next@ --- *
154 *
155 * Arguments: @primeiter *pi@ = pointer to an iterator structure
156 * @mp *d@ = fake destination
157 *
158 * Returns: The next prime number from the iterator.
159 *
160 * Use: Returns a new prime number.
161 */
162
163mp *primeiter_next(primeiter *pi, mp *d)
164{
165 mp *p;
166
167 switch (pi->mode) {
168 case PIMODE_PTAB:
169 pi->w = primetab[pi->i++];
170 if (pi->i >= NPRIME) {
171 wheelsync(pi, pi->p);
172 pi->mode = PIMODE_WHEEL;
173 }
174 p = MP_COPY(pi->p);
175 MP_SPLIT(p);
176 break;
177 case PIMODE_STALL:
178 pi->mode = PIMODE_WHEEL;
179 goto loop;
180 case PIMODE_WHEEL:
181 do {
182 MP_DEST(pi->p, MP_LEN(pi->p) + 1, pi->p->f);
183 MPX_UADDN(pi->p->v, pi->p->vl, wheel[pi->i++]);
184 MP_SHRINK(pi->p);
185 if (pi->i >= WHEELN) pi->i = 0;
186 loop:;
187 } while (!pgen_primep(pi->p, pi->r));
188 p = MP_COPY(pi->p);
189 break;
190 default:
191 abort();
192 }
193 if (d) MP_DROP(d);
194 return (p);
195}
196
197/*----- Test rig ----------------------------------------------------------*/
198
199#ifdef TEST_RIG
200
201#include <mLib/macros.h>
202#include <mLib/testrig.h>
203
204static int test(dstr *v)
205{
206 mp *start = *(mp **)v[0].buf;
207 mp *pp[5], *ret[5];
208 int i;
209 primeiter pi;
210 int ok = 1;
211
212 for (i = 0; i < N(pp); i++)
213 pp[i] = *(mp **)v[i + 1].buf;
214 primeiter_create(&pi, start);
215 for (i = 0; i < N(pp); i++) {
216 ret[i] = primeiter_next(&pi, MP_NEW);
217 if (!MP_EQ(ret[i], pp[i])) ok = 0;
218 }
219 if (!ok) {
220 fprintf(stderr, "\n*** primeiter test failure:\n*** start = ");
221 mp_writefile(start, stderr, 10);
222 for (i = 0; i < N(pp); i++) {
223 fprintf(stderr, "\n*** p[%d] = ", i);
224 mp_writefile(ret[i], stderr, 10);
225 fprintf(stderr, " %s ", MP_EQ(ret[i], pp[i]) ? "==" : "!=");
226 mp_writefile(pp[i], stderr, 10);
227 }
228 fputc('\n', stderr);
229 }
230 for (i = 0; i < N(pp); i++) {
231 MP_DROP(pp[i]);
232 MP_DROP(ret[i]);
233 }
234 primeiter_destroy(&pi);
235 MP_DROP(start);
236 assert(mparena_count(MPARENA_GLOBAL) == 0);
237 return (ok);
238}
239
240static test_chunk tests[] = {
241 { "primeiter", test,
242 { &type_mp, &type_mp, &type_mp, &type_mp, &type_mp, &type_mp, } },
243 { 0 }
244};
245
246int main(int argc, char *argv[])
247{
0f00dc4c 248 test_run(argc, argv, tests, SRCDIR "/t/pgen");
1d6d3b01
MW
249 return (0);
250}
251
252#endif
253
254/*----- That's all, folks -------------------------------------------------*/