3 * Prime generation glue
5 * (c) 1999 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 ------------------------------------------------------*/
43 /*----- Standard prime filter ---------------------------------------------*/
45 /* --- @pgen_filter@ --- */
47 int pgen_filter(int rq, pgen_event *ev, void *p)
49 pgen_filterctx *f = p;
54 rc = pfilt_create(&f->f, ev->m);
68 if (rc == PGEN_FAIL && !((f->step | f->f.m->v[0]) & 1))
69 rc = pfilt_step(&f->f, 1);
70 while (rc == PGEN_FAIL)
71 rc = pfilt_step(&f->f, f->step);
72 ev->m = MP_COPY(f->f.m);
76 /* --- @pgen_jump@ --- *
78 * Similar to the standard @pgen_filter@, but jumps in large steps rather
82 int pgen_jump(int rq, pgen_event *ev, void *p)
90 mp_gcd(&g, 0, 0, ev->m, f->j->m);
91 if (MP_CMP(g, >, MP_ONE)) {
96 rc = pfilt_create(&f->f, ev->m);
101 rc = pfilt_jump(&f->f, f->j);
104 pfilt_destroy(&f->f);
108 while (rc == PGEN_FAIL)
109 rc = pfilt_jump(&f->f, f->j);
110 ev->m = MP_COPY(f->f.m);
114 /*----- Standard prime test -----------------------------------------------*/
116 /* --- @pgen_test@ --- */
118 int pgen_test(int rq, pgen_event *ev, void *p)
125 rabin_create(r, ev->m);
130 rc = rabin_rtest(r, MP_TWO);
132 mp *a = mprand_range(MP_NEW, ev->m, ev->r, 0);
133 rc = rabin_rtest(r, a);
146 /*----- The main driver ---------------------------------------------------*/
150 * Arguments: @const char *name@ = name of the value being searched for
151 * @mp *d@ = destination for the result integer
152 * @mp *m@ = start value to pass to stepper
153 * @pgen_proc *event@ = event handler function
154 * @void *ectx@ = context argument for event andler
155 * @unsigned steps@ = number of steps to take in search
156 * @pgen_proc *step@ = stepper function to use
157 * @void *sctx@ = context argument for stepper
158 * @unsigned tests@ = number of tests to make
159 * @pgen_proc *test@ = tester function to use
160 * @void *tctx@ = context argument for tester
162 * Returns: Pointer to final result, or null.
164 * Use: A generalized prime-number search skeleton. Yes, that's a
165 * scary number of arguments.
168 mp *pgen(const char *name, mp *d, mp *m, pgen_proc *event, void *ectx,
169 unsigned steps, pgen_proc *step, void *sctx,
170 unsigned tests, pgen_proc *test, void *tctx)
178 enum { P_STEP, P_TEST };
180 /* --- Set up the initial event block --- */
189 ev.r = fibrand_create(0);
191 /* --- Tell the event handler we're under way --- */
193 if (event && event(PGEN_BEGIN, &ev, ectx) == PGEN_ABORT) {
194 ev.r->ops->destroy(ev.r);
198 /* --- Set up for the initial call --- */
200 proc = step; ctx = sctx; p = P_STEP; rq = PGEN_BEGIN;
202 /* --- Enter the great maelstrom of state transitions --- */
211 #define A_ENDSTEP 16u
214 /* --- Call the procedure and decide what to do next --- */
216 rc = proc(rq, &ev, ctx);
223 proc = test; ctx = tctx; p = P_TEST;
228 act |= A_TEST | A_EVENT;
232 proc = test; ctx = tctx; p = P_TEST;
239 act |= A_ENDTEST | A_EVENT;
240 proc = step; ctx = sctx; p = P_STEP;
245 act |= A_EVENT | A_DONE | A_ENDSTEP;
250 act |= A_EVENT | A_DONE;
251 if (p == P_TEST || rq == PGEN_TRY)
253 if (p == P_TEST && rq != PGEN_BEGIN)
257 assert(((void)"Invalid response from function", 0));
261 /* --- If decrementing counters is requested, do that --- */
263 if ((act & A_STEP) && steps) {
265 if (ev.steps == steps) {
266 act |= A_EVENT | A_ENDSTEP | A_DONE;
272 if ((act & A_TEST) && tests) {
274 if (ev.tests == tests) {
275 act |= A_ENDTEST | A_ENDSTEP | A_DONE;
280 /* --- Report an event if so directed --- */
282 if ((act & A_EVENT) && event && event(rc, &ev, ectx) == PGEN_ABORT) {
284 if (!(act & A_DONE)) {
285 act |= A_ENDSTEP | A_DONE;
286 if (p == P_TEST && rq != PGEN_BEGIN)
291 /* --- Close down tester and stepper functions --- */
294 test(PGEN_DONE, &ev, tctx);
296 step(PGEN_DONE, &ev, sctx);
298 /* --- Stop the entire test if necessary --- */
304 /* --- Tidy up and return --- */
306 if (rc == PGEN_ABORT) {
310 ev.r->ops->destroy(ev.r);
316 /* --- @pgen_primep@ --- *
318 * Arguments: @mp *p@ = a number to check
319 * @grand *gr@ = a random number source
321 * Returns: Nonzero if @p@ is really prime.
323 * Use: Checks the primality of @p@. If @p@ is prime, then this
324 * function returns nonzero; if @p@ is really composite then it
325 * %%\emph{probably}%% returns zero, but might not.
327 * Currently, this function uses the Baillie--PSW test, which
328 * combines a single Miller--Rabin test with witness 2 with a
329 * single Frobenius test with parameters chosen using
330 * Selfridge's `Method A'. No composites are known which pass
331 * this test, though it's conjectured that infinitely many
335 int pgen_primep(mp *p, grand *gr)
340 if (MP_NEGP(p)) return (0);
341 switch (pfilt_smallfactor(p)) {
342 case PGEN_DONE: return (1);
343 case PGEN_FAIL: return (0);
345 rabin_create(&r, p); rc = rabin_test(&r, MP_TWO); rabin_destroy(&r);
346 if (rc == PGEN_FAIL) return (0);
347 rc = pgen_granfrob(p, 0, 0); if (rc == PGEN_FAIL) return (0);
351 /*----- Test rig ----------------------------------------------------------*/
355 #include <mLib/testrig.h>
357 static int t_primep(dstr *v)
359 mp *m = *(mp **)v[0].buf;
360 int e = *(int *)v[1].buf;
365 rng = fibrand_create(0);
366 r = pgen_primep(m, rng);
369 fputs("\n*** primep failed", stderr);
370 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
371 fprintf(stderr, "\nexpected %d", e);
372 fprintf(stderr, "\nreported %d", r);
378 assert(mparena_count(MPARENA_GLOBAL) == 0);
382 static int verify(dstr *v)
384 mp *m = *(mp **)v[0].buf;
385 mp *q = *(mp **)v[1].buf;
393 p = pgen("p", MP_NEW, m, pgen_evspin, 0, 0, pgen_filter, &pf,
394 rabin_iters(mp_bits(m)), pgen_test, &r);
395 if (!p || !MP_EQ(p, q)) {
396 fputs("\n*** pgen failed", stderr);
397 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
398 fputs("\np = ", stderr); mp_writefile(p, stderr, 10);
399 fputs("\nq = ", stderr); mp_writefile(q, stderr, 10);
407 assert(mparena_count(MPARENA_GLOBAL) == 0);
411 static test_chunk tests[] = {
412 { "pgen", verify, { &type_mp, &type_mp, 0 } },
413 { "primep", t_primep, { &type_mp, &type_int, 0 } },
417 int main(int argc, char *argv[])
420 test_run(argc, argv, tests, SRCDIR "/t/pgen");
425 /*----- That's all, folks -------------------------------------------------*/