3 * Iterate over small primes efficiently
5 * (c) 2007 Straylight/Edgeware
8 /*----- Licensing notice --------------------------------------------------*
10 * This file is part of Catacomb.
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.
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.
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,
28 /*----- Header files ------------------------------------------------------*/
33 #include "primeiter.h"
37 /*----- Theory ------------------------------------------------------------*
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.
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.
52 /*----- Main code ---------------------------------------------------------*/
54 /* --- @wheelsync@ --- *
56 * Arguments: @primeiter *pi@ = iterator to synchronize
57 * @mp *where@ = value to synchronize
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.
65 static void wheelsync(primeiter
*pi
, mp
*where
)
74 mp_build(&t
, &w
, &w
+ 1);
75 mp_div(0, &r
, where
, &t
);
76 rr
= MP_ZEROP(r
) ?
0 : r
->v
[0];
78 for (i
= 0, n
= 1; rr
> n
; n
+= wheel
[i
], i
++);
80 pi
->p
= mp_add(MP_NEW
, where
, &t
);
82 pi
->r
= fibrand_create(0);
86 /* --- @primeiter_create@ --- *
88 * Arguments: @primeiter *pi@ = pointer to an iterator structure
89 * @mp *start@ = where to start
93 * Use: Initializes a prime iterator. The first output will be the
94 * smallest prime not less than @start@.
97 void primeiter_create(primeiter
*pi
, mp
*start
)
102 if (!start
|| MP_CMP(start
, <=, MP_TWO
))
105 if (MP_LEN(start
) <= 1) {
112 if (primetab
[m
] == n
) break;
113 else if (m
== l
) { m
++; break; }
114 else if (primetab
[m
] < n
) l
= m
;
118 pi
->mode
= PIMODE_PTAB
;
119 mp_build(&pi
->pp
, &pi
->w
, &pi
->w
+ 1);
125 wheelsync(pi
, start
);
126 pi
->mode
= PIMODE_STALL
;
129 /* --- @primeiter_destroy@ --- *
131 * Arguments: @primeiter *pi@ = pointer to iterator structure
135 * Use: Frees up an iterator structure when it's no longer wanted.
138 void primeiter_destroy(primeiter
*pi
)
153 /* --- @primeiter_next@ --- *
155 * Arguments: @primeiter *pi@ = pointer to an iterator structure
156 * @mp *d@ = fake destination
158 * Returns: The next prime number from the iterator.
160 * Use: Returns a new prime number.
163 mp
*primeiter_next(primeiter
*pi
, mp
*d
)
169 pi
->w
= primetab
[pi
->i
++];
170 if (pi
->i
>= NPRIME
) {
171 wheelsync(pi
, pi
->p
);
172 pi
->mode
= PIMODE_WHEEL
;
178 pi
->mode
= PIMODE_WHEEL
;
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
++]);
185 if (pi
->i
>= WHEELN
) pi
->i
= 0;
187 } while (!pgen_primep(pi
->p
, pi
->r
));
197 /*----- Test rig ----------------------------------------------------------*/
201 #include <mLib/macros.h>
202 #include <mLib/testrig.h>
204 static int test(dstr
*v
)
206 mp
*start
= *(mp
**)v
[0].buf
;
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;
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);
230 for (i
= 0; i
< N(pp
); i
++) {
234 primeiter_destroy(&pi
);
236 assert(mparena_count(MPARENA_GLOBAL
) == 0);
240 static test_chunk tests
[] = {
242 { &type_mp
, &type_mp
, &type_mp
, &type_mp
, &type_mp
, &type_mp
, } },
246 int main(int argc
, char *argv
[])
248 test_run(argc
, argv
, tests
, SRCDIR
"/t/pgen");
254 /*----- That's all, folks -------------------------------------------------*/