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
)
126 rabin_create(r
, ev
->m
);
130 a
= mprand_range(a
, ev
->m
, ev
->r
, 0);
131 rc
= rabin_rtest(r
, a
);
143 /* --- @pgen_bailliepswtest@ --- */
145 int pgen_bailliepswtest(int rq
, pgen_event
*ev
, void *p
)
152 if (ev
->tests
!= 2) rc
= PGEN_ABORT
;
163 rabin_create(&r
, ev
->m
);
164 rc
= rabin_test(&r
, MP_TWO
);
168 rc
= pgen_granfrob(ev
->m
, 0, 0);
184 /*----- The main driver ---------------------------------------------------*/
188 * Arguments: @const char *name@ = name of the value being searched for
189 * @mp *d@ = destination for the result integer
190 * @mp *m@ = start value to pass to stepper
191 * @pgen_proc *event@ = event handler function
192 * @void *ectx@ = context argument for event andler
193 * @unsigned steps@ = number of steps to take in search
194 * @pgen_proc *step@ = stepper function to use
195 * @void *sctx@ = context argument for stepper
196 * @unsigned tests@ = number of tests to make
197 * @pgen_proc *test@ = tester function to use
198 * @void *tctx@ = context argument for tester
200 * Returns: Pointer to final result, or null.
202 * Use: A generalized prime-number search skeleton. Yes, that's a
203 * scary number of arguments.
206 mp
*pgen(const char *name
, mp
*d
, mp
*m
, pgen_proc
*event
, void *ectx
,
207 unsigned steps
, pgen_proc
*step
, void *sctx
,
208 unsigned tests
, pgen_proc
*test
, void *tctx
)
216 enum { P_STEP
, P_TEST
};
218 /* --- Set up the initial event block --- */
227 ev
.r
= fibrand_create(0);
229 /* --- Tell the event handler we're under way --- */
231 if (event
&& event(PGEN_BEGIN
, &ev
, ectx
) == PGEN_ABORT
) {
232 ev
.r
->ops
->destroy(ev
.r
);
236 /* --- Set up for the initial call --- */
238 proc
= step
; ctx
= sctx
; p
= P_STEP
; rq
= PGEN_BEGIN
;
240 /* --- Enter the great maelstrom of state transitions --- */
249 #define A_ENDSTEP 16u
252 /* --- Call the procedure and decide what to do next --- */
254 rc
= proc(rq
, &ev
, ctx
);
261 proc
= test
; ctx
= tctx
; p
= P_TEST
;
266 act
|= A_TEST
| A_EVENT
;
270 proc
= test
; ctx
= tctx
; p
= P_TEST
;
277 act
|= A_ENDTEST
| A_EVENT
;
278 proc
= step
; ctx
= sctx
; p
= P_STEP
;
283 act
|= A_EVENT
| A_DONE
| A_ENDSTEP
;
288 act
|= A_EVENT
| A_DONE
;
289 if (p
== P_TEST
|| rq
== PGEN_TRY
)
291 if (p
== P_TEST
&& rq
!= PGEN_BEGIN
)
295 assert(((void)"Invalid response from function", 0));
299 /* --- If decrementing counters is requested, do that --- */
301 if ((act
& A_STEP
) && steps
) {
304 act
|= A_EVENT
| A_ENDSTEP
| A_DONE
;
310 if ((act
& A_TEST
) && tests
) {
313 act
|= A_ENDTEST
| A_ENDSTEP
| A_DONE
;
318 /* --- Report an event if so directed --- */
320 if ((act
& A_EVENT
) && event
&& event(rc
, &ev
, ectx
) == PGEN_ABORT
) {
322 if (!(act
& A_DONE
)) {
323 act
|= A_ENDSTEP
| A_DONE
;
324 if (p
== P_TEST
&& rq
!= PGEN_BEGIN
)
329 /* --- Close down tester and stepper functions --- */
332 test(PGEN_DONE
, &ev
, tctx
);
334 step(PGEN_DONE
, &ev
, sctx
);
336 /* --- Stop the entire test if necessary --- */
342 /* --- Tidy up and return --- */
344 if (rc
== PGEN_ABORT
) {
348 ev
.r
->ops
->destroy(ev
.r
);
354 /* --- @pgen_primep@ --- *
356 * Arguments: @mp *p@ = a number to check
357 * @grand *gr@ = a random number source
359 * Returns: Nonzero if @p@ is really prime.
361 * Use: Checks the primality of @p@. If @p@ is prime, then this
362 * function returns nonzero; if @p@ is really composite then it
363 * %%\emph{probably}%% returns zero, but might not.
365 * Currently, this function uses the Baillie--PSW test, which
366 * combines a single Miller--Rabin test with witness 2 with a
367 * single Frobenius test with parameters chosen using
368 * Selfridge's `Method A'. No composites are known which pass
369 * this test, though it's conjectured that infinitely many
373 int pgen_primep(mp
*p
, grand
*gr
)
378 if (MP_NEGP(p
)) return (0);
379 switch (pfilt_smallfactor(p
)) {
380 case PGEN_DONE
: return (1);
381 case PGEN_FAIL
: return (0);
383 rabin_create(&r
, p
); rc
= rabin_test(&r
, MP_TWO
); rabin_destroy(&r
);
384 if (rc
== PGEN_FAIL
) return (0);
385 rc
= pgen_granfrob(p
, 0, 0); if (rc
== PGEN_FAIL
) return (0);
389 /*----- Test rig ----------------------------------------------------------*/
393 #include <mLib/testrig.h>
395 static int t_primep(dstr
*v
)
397 mp
*m
= *(mp
**)v
[0].buf
;
398 int e
= *(int *)v
[1].buf
;
403 rng
= fibrand_create(0);
404 r
= pgen_primep(m
, rng
);
407 fputs("\n*** primep failed", stderr
);
408 fputs("\nm = ", stderr
); mp_writefile(m
, stderr
, 10);
409 fprintf(stderr
, "\nexpected %d", e
);
410 fprintf(stderr
, "\nreported %d", r
);
416 assert(mparena_count(MPARENA_GLOBAL
) == 0);
420 static int verify(dstr
*v
)
422 mp
*m
= *(mp
**)v
[0].buf
;
423 mp
*q
= *(mp
**)v
[1].buf
;
431 p
= pgen("p", MP_NEW
, m
, pgen_evspin
, 0, 0, pgen_filter
, &pf
,
432 rabin_iters(mp_bits(m
)), pgen_test
, &r
);
433 if (!p
|| !MP_EQ(p
, q
)) {
434 fputs("\n*** pgen failed", stderr
);
435 fputs("\nm = ", stderr
); mp_writefile(m
, stderr
, 10);
436 fputs("\np = ", stderr
); mp_writefile(p
, stderr
, 10);
437 fputs("\nq = ", stderr
); mp_writefile(q
, stderr
, 10);
445 assert(mparena_count(MPARENA_GLOBAL
) == 0);
449 static test_chunk tests
[] = {
450 { "pgen", verify
, { &type_mp
, &type_mp
, 0 } },
451 { "primep", t_primep
, { &type_mp
, &type_int
, 0 } },
455 int main(int argc
, char *argv
[])
458 test_run(argc
, argv
, tests
, SRCDIR
"/t/pgen");
463 /*----- That's all, folks -------------------------------------------------*/