X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/d94f85acc2dce2b32cd7d2870a749d20dd263946..5278d9afdf1aff9fd6f64073ea42395d756ee58c:/pgen.c diff --git a/pgen.c b/pgen.c index e9a5e75..1a2a8ac 100644 --- a/pgen.c +++ b/pgen.c @@ -1,13 +1,13 @@ /* -*-c-*- * - * $Id: pgen.c,v 1.3 1999/12/10 23:28:35 mdw Exp $ + * $Id$ * - * Finding and testing prime numbers + * Prime generation glue * * (c) 1999 Straylight/Edgeware */ -/*----- Licensing notice --------------------------------------------------* +/*----- Licensing notice --------------------------------------------------* * * This file is part of Catacomb. * @@ -15,324 +15,399 @@ * it under the terms of the GNU Library General Public License as * published by the Free Software Foundation; either version 2 of the * License, or (at your option) any later version. - * + * * Catacomb is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Library General Public License for more details. - * + * * You should have received a copy of the GNU Library General Public * License along with Catacomb; if not, write to the Free * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, * MA 02111-1307, USA. */ -/*----- Revision history --------------------------------------------------* - * - * $Log: pgen.c,v $ - * Revision 1.3 1999/12/10 23:28:35 mdw - * Track suggested destination changes. - * - * Revision 1.2 1999/11/20 22:23:05 mdw - * Add multiply-and-add function for Diffie-Hellman safe prime generation. - * - * Revision 1.1 1999/11/19 13:17:57 mdw - * Prime number generator and tester. - * - */ - /*----- Header files ------------------------------------------------------*/ +#include +#include +#include +#include + +#include "fibrand.h" +#include "grand.h" #include "mp.h" -#include "mpmont.h" +#include "mprand.h" #include "pgen.h" -#include "ptab.h" +#include "pfilt.h" #include "rabin.h" -/*----- Main code ---------------------------------------------------------*/ +/*----- Standard prime filter ---------------------------------------------*/ -/* --- @pgen_create@ --- * - * - * Arguments: @pgen *p@ = pointer to prime generation context - * @mp *m@ = pointer to initial number to test - * - * Returns: One of the @PGEN@ constants above. - * - * Use: Tests an initial number for primality by computing its - * residue modulo various small prime numbers. This is fairly - * quick, but not particularly certain. If a @PGEN_MAYBE@ - * result is returned, perform Rabin-Miller tests to confirm. - */ +/* --- @pgen_filter@ --- */ -int pgen_create(pgen *p, mp *m) +int pgen_filter(int rq, pgen_event *ev, void *p) { - int rc = PGEN_MAYBE; - int i; - mp *r = MP_NEW; - mpw qw; - mp q; - - /* --- Take a copy of the number --- */ - - mp_shrink(m); - p->m = MP_COPY(m); - - /* --- Fill in the residues --- */ - - mp_build(&q, &qw, &qw + 1); - for (i = 0; i < NPRIME; i++) { - qw = ptab[i]; - mp_div(0, &r, m, &q); - p->r[i] = r->v[0]; - if (!p->r[i] && rc == PGEN_MAYBE) { - if (MP_LEN(m) == 1 && m->v[0] == ptab[i]) - rc = PGEN_PRIME; - else - rc = PGEN_COMPOSITE; - } - } + pgen_filterctx *f = p; + int rc = PGEN_FAIL; - /* --- Done --- */ + switch (rq) { + case PGEN_BEGIN: + rc = pfilt_create(&f->f, ev->m); + mp_drop(ev->m); + break; + case PGEN_TRY: + mp_drop(ev->m); + break; + case PGEN_DONE: + pfilt_destroy(&f->f); + return (PGEN_DONE); + default: + rc = PGEN_ABORT; + break; + } - mp_drop(r); + if (rc == PGEN_FAIL && !((f->step | f->f.m->v[0]) & 1)) + rc = pfilt_step(&f->f, 1); + while (rc == PGEN_FAIL) + rc = pfilt_step(&f->f, f->step); + ev->m = MP_COPY(f->f.m); return (rc); } -/* --- @pgen_destroy@ --- * - * - * Arguments: @pgen *p@ = pointer to prime generation context - * - * Returns: --- +/* --- @pgen_jump@ --- * * - * Use: Discards a context and all the resources it holds. + * Similar to the standard @pgen_filter@, but jumps in large steps rather + * than small ones. */ -void pgen_destroy(pgen *p) +int pgen_jump(int rq, pgen_event *ev, void *p) { - mp_drop(p->m); -} - -/* --- @pgen_step@ --- * - * - * Arguments: @pgen *p@ = pointer to prime generation context - * @mpw step@ = how much to step the number - * - * Returns: One of the @PGEN@ constants above. - * - * Use: Steps a number by a small amount. Stepping is much faster - * than initializing with a new number. The test performed is - * the same simple one used by @ptab_create@, so @PGEN_MAYBE@ - * results should be followed up by a Rabin-Miller test. - */ + pgen_jumpctx *f = p; + int rc = PGEN_ABORT; + + switch (rq) { + case PGEN_BEGIN: { + mp *g = MP_NEW; + mp_gcd(&g, 0, 0, ev->m, f->j->m); + if (MP_CMP(g, >, MP_ONE)) { + mp_drop(g); + return (PGEN_ABORT); + } + mp_drop(g); + rc = pfilt_create(&f->f, ev->m); + mp_drop(ev->m); + } break; + case PGEN_TRY: + mp_drop(ev->m); + rc = pfilt_jump(&f->f, f->j); + break; + case PGEN_DONE: + pfilt_destroy(&f->f); + return (PGEN_DONE); + } -int pgen_step(pgen *p, mpw step) -{ - int rc = PGEN_MAYBE; - int i; + while (rc == PGEN_FAIL) + rc = pfilt_jump(&f->f, f->j); + ev->m = MP_COPY(f->f.m); + return (rc); +} - /* --- Add the step on to the number --- */ +/*----- Standard prime test -----------------------------------------------*/ - p->m = mp_split(p->m); - mp_ensure(p->m, MP_LEN(p->m) + 1); - mpx_uaddn(p->m->v, p->m->vl, step); - mp_shrink(p->m); +/* --- @pgen_test@ --- */ - /* --- Update the residue table --- */ +int pgen_test(int rq, pgen_event *ev, void *p) +{ + rabin *r = p; + int rc = PGEN_ABORT; - for (i = 0; i < NPRIME; i++) { - p->r[i] = (p->r[i] + step) % ptab[i]; - if (!p->r[i] && rc == PGEN_MAYBE) { - if (MP_LEN(p->m) == 1 && p->m->v[0] == ptab[i]) - rc = PGEN_PRIME; - else - rc = PGEN_COMPOSITE; - } + switch (rq) { + case PGEN_BEGIN: + rabin_create(r, ev->m); + rc = PGEN_TRY; + break; + case PGEN_TRY: + if (!ev->tests) + rc = rabin_rtest(r, MP_TWO); + else { + mp *a = mprand_range(MP_NEW, ev->m, ev->r, 0); + rc = rabin_rtest(r, a); + mp_drop(a); + } + break; + case PGEN_DONE: + rabin_destroy(r); + rc = PGEN_DONE; + break; } - /* --- Small numbers must be prime --- */ - - if (rc == PGEN_MAYBE && MP_LEN(p->m) == 1 && - p->m->v[0] < MAXPRIME * MAXPRIME) - rc = PGEN_PRIME; - - /* --- Done --- */ - return (rc); } -/* --- @pgen_muladd@ --- * - * - * Arguments: @pgen *p@ = destination prime generation context - * @const pgen *q@ = source prime generation context - * @mpw m@ = number to multiply by - * @mpw a@ = number to add +/*----- The main driver ---------------------------------------------------*/ + +/* --- @pgen@ --- * * - * Returns: One of the @PGEN@ constants above. + * Arguments: @const char *name@ = name of the value being searched for + * @mp *d@ = destination for the result integer + * @mp *m@ = start value to pass to stepper + * @pgen_proc *event@ = event handler function + * @void *ectx@ = context argument for event andler + * @unsigned steps@ = number of steps to take in search + * @pgen_proc *step@ = stepper function to use + * @void *sctx@ = context argument for stepper + * @unsigned tests@ = number of tests to make + * @pgen_proc *test@ = tester function to use + * @void *tctx@ = context argument for tester * - * Use: Multiplies the number in a prime generation context by a - * small value and then adds a small value. The destination - * should either be uninitialized or the same as the source. + * Returns: Pointer to final result, or null. * - * Common things to do include multiplying by 2 and adding 0 to - * turn a prime into a jump for finding other primes with @q@ as - * a factor of @p - 1@, or multiplying by 2 and adding 1. + * Use: A generalized prime-number search skeleton. Yes, that's a + * scary number of arguments. */ -int pgen_muladd(pgen *p, const pgen *q, mpw m, mpw a) +mp *pgen(const char *name, mp *d, mp *m, pgen_proc *event, void *ectx, + unsigned steps, pgen_proc *step, void *sctx, + unsigned tests, pgen_proc *test, void *tctx) { - int rc = PGEN_MAYBE; - int i; - - /* --- Multiply the big number --- */ - - { - mp *d = mp_create(MP_LEN(q->m) + 2); - mpx_umuln(d->v, d->vl, q->m->v, q->m->vl, m); - mpx_uaddn(d->v, d->vl, a); - d->f = q->m->f; - if (p == q) - mp_drop(p->m); - mp_shrink(d); - p->m = d; + pgen_event ev; + int rq, rc; + pgen_proc *proc; + void *ctx; + int p; + + enum { P_STEP, P_TEST }; + + /* --- Set up the initial event block --- */ + + ev.name = name; + if (m) + ev.m = MP_COPY(m); + else + ev.m = 0; + ev.steps = 0; + ev.tests = 0; + ev.r = fibrand_create(0); + + /* --- Tell the event handler we're under way --- */ + + if (event && event(PGEN_BEGIN, &ev, ectx) == PGEN_ABORT) { + ev.r->ops->destroy(ev.r); + return (0); } - /* --- Gallivant through the residue table --- */ - - for (i = 0; i < NPRIME; i++) { - p->r[i] = (q->r[i] * m + a) % ptab[i]; - if (!p->r[i] && rc == PGEN_MAYBE) { - if (MP_LEN(p->m) == 1 && p->m->v[0] == ptab[i]) - rc = PGEN_PRIME; - else - rc = PGEN_COMPOSITE; + /* --- Set up for the initial call --- */ + + proc = step; ctx = sctx; p = P_STEP; rq = PGEN_BEGIN; + + /* --- Enter the great maelstrom of state transitions --- */ + + for (;;) { + unsigned act = 0; + +#define A_STEP 1u +#define A_TEST 2u +#define A_EVENT 4u +#define A_ENDTEST 8u +#define A_ENDSTEP 16u +#define A_DONE 32u + + /* --- Call the procedure and decide what to do next --- */ + + rc = proc(rq, &ev, ctx); + switch (rc) { + case PGEN_TRY: + if (p == P_TEST) + rq = PGEN_TRY; + else { + act |= A_EVENT; + proc = test; ctx = tctx; p = P_TEST; + rq = PGEN_BEGIN; + } + break; + case PGEN_PASS: + act |= A_TEST | A_EVENT; + if (p == P_TEST) + rq = PGEN_TRY; + else { + proc = test; ctx = tctx; p = P_TEST; + rq = PGEN_BEGIN; + } + break; + case PGEN_FAIL: + act |= A_STEP; + if (p == P_TEST) { + act |= A_ENDTEST | A_EVENT; + proc = step; ctx = sctx; p = P_STEP; + } + rq = PGEN_TRY; + break; + case PGEN_DONE: + act |= A_EVENT | A_DONE | A_ENDSTEP; + if (p == P_TEST) + act |= A_ENDTEST; + break; + case PGEN_ABORT: + act |= A_EVENT | A_DONE; + if (p == P_TEST || rq == PGEN_TRY) + act |= A_ENDSTEP; + if (p == P_TEST && rq != PGEN_BEGIN) + act |= A_ENDTEST; + break; + default: + assert(((void)"Invalid response from function", 0)); + break; } - } - /* --- Small numbers must be prime --- */ + /* --- If decrementing counters is requested, do that --- */ - if (rc == PGEN_MAYBE && MP_LEN(p->m) == 1 && - p->m->v[0] < MAXPRIME * MAXPRIME) - rc = PGEN_PRIME; + if ((act & A_STEP) && steps) { + ev.steps++; + if (ev.steps == steps) { + act |= A_EVENT | A_ENDSTEP | A_DONE; + rc = PGEN_ABORT; + } + ev.tests = 0; + } - /* --- Finished --- */ + if ((act & A_TEST) && tests) { + ev.tests++; + if (ev.tests == tests) { + act |= A_ENDTEST | A_ENDSTEP | A_DONE; + rc = PGEN_DONE; + } + } - return (rc); -} + /* --- Report an event if so directed --- */ + if ((act & A_EVENT) && event && event(rc, &ev, ectx) == PGEN_ABORT) { + rc = PGEN_ABORT; + if (!(act & A_DONE)) { + act |= A_ENDSTEP | A_DONE; + if (p == P_TEST) + act |= A_ENDTEST; + } + } -/* --- @pgen_jump@ --- * - * - * Arguments: @pgen *p@ = pointer to prime generation context - * @pgen *j@ = pointer to another generation context - * - * Returns: One of the @PGEN@ constants above. - * - * Use: Steps a number by a large amount. Even so, jumping is much - * faster than initializing a new number. The test peformed is - * the same simple one used by @ptab_create@, so @PGEN_MAYBE@ - * results should be followed up by a Rabin-Miller test. - * - * Note that the number stored in the @j@ context is probably - * better off being even than prime. The important thing is - * that all of the residues for the number have already been - * computed. - */ + /* --- Close down tester and stepper functions --- */ -int pgen_jump(pgen *p, pgen *j) -{ - int rc = PGEN_MAYBE; - int i; + if (act & A_ENDTEST) + test(PGEN_DONE, &ev, tctx); + if (act & A_ENDSTEP) + step(PGEN_DONE, &ev, sctx); - /* --- Add the step on --- */ + /* --- Stop the entire test if necessary --- */ - p->m = mp_add(p->m, p->m, j->m); + if (act & A_DONE) + break; + } - /* --- Update the residue table --- */ + /* --- Tidy up and return --- */ - for (i = 0; i < NPRIME; i++) { - p->r[i] = p->r[i] + j->r[i]; - if (p->r[i] > ptab[i]) - p->r[i] -= ptab[i]; - if (!p->r[i] && rc == PGEN_MAYBE) { - if (MP_LEN(p->m) == 1 && p->m->v[0] == ptab[i]) - rc = PGEN_PRIME; - else - rc = PGEN_COMPOSITE; - } + if (rc == PGEN_ABORT) { + mp_drop(ev.m); + ev.m = 0; } + ev.r->ops->destroy(ev.r); + mp_drop(d); - /* --- Small numbers must be prime --- */ + return (ev.m); +} - if (rc == PGEN_MAYBE && MP_LEN(p->m) == 1 && - p->m->v[0] < MAXPRIME * MAXPRIME) - rc = PGEN_PRIME; +/* --- @pgen_primep@ --- * + * + * Arguments: @mp *p@ = a number to check + * @grand *gr@ = a random number source + * + * Returns: Nonzero if @p@ is really prime. + */ - /* --- Done --- */ +int pgen_primep(mp *p, grand *gr) +{ + int i; + rabin r; + mp *x = MP_NEW; - return (rc); + if (MP_NEGP(p)) return (0); + switch (pfilt_smallfactor(p)) { + case PGEN_DONE: return (1); + case PGEN_FAIL: return (0); + } + rabin_create(&r, p); + for (i = 32; i; i--) { + x = mprand_range(x, p, gr, 0); + if (rabin_rtest(&r, x) == PGEN_FAIL) + break; + } + MP_DROP(x); + rabin_destroy(&r); + return (!i); } -/*----- Test code ---------------------------------------------------------*/ +/*----- Test rig ----------------------------------------------------------*/ #ifdef TEST_RIG #include -static int verify(dstr *v) +static int t_primep(dstr *v) { mp *m = *(mp **)v[0].buf; - mp *r = *(mp **)v[1].buf; - pgen p; - mpw gw; - mp g; - int rc; - static char baton[] = "/-\\|"; - char *b = baton; + int e = *(int *)v[1].buf; + int r; + grand *rng; int ok = 1; - mp_build(&g, &gw, &gw + 1); - rc = pgen_create(&p, m); - for (;;) { - if (rc == PGEN_PRIME) - break; - if (rc == PGEN_MAYBE) { - int i; - rabin r; - rabin_create(&r, p.m); - for (i = 0; i < 5; i++) { - gw = rand(); - if (rabin_test(&r, &g) == PGEN_COMPOSITE) - break; - } - rabin_destroy(&r); - if (i == 5) - break; - putc(*b++, stderr); - putc('\b', stderr); - if (!*b) - b = baton; - } - rc = pgen_step(&p, 2); + rng = fibrand_create(0); + r = pgen_primep(m, rng); + GR_DESTROY(rng); + if (e != r) { + fputs("\n*** primep failed", stderr); + fputs("\nm = ", stderr); mp_writefile(m, stderr, 10); + fprintf(stderr, "\nexpected %d", e); + fprintf(stderr, "\nreported %d", r); + fputc('\n', stderr); + ok = 0; } - if (MP_CMP(p.m, !=, r)) { - fputs("\n*** failed.", stderr); + mp_drop(m); + assert(mparena_count(MPARENA_GLOBAL) == 0); + return (ok); +} + +static int verify(dstr *v) +{ + mp *m = *(mp **)v[0].buf; + mp *q = *(mp **)v[1].buf; + mp *p; + int ok = 1; + + pgen_filterctx pf; + rabin r; + + pf.step = 2; + p = pgen("p", MP_NEW, m, pgen_evspin, 0, 0, pgen_filter, &pf, + rabin_iters(mp_bits(m)), pgen_test, &r); + if (!p || !MP_EQ(p, q)) { + fputs("\n*** pgen failed", stderr); fputs("\nm = ", stderr); mp_writefile(m, stderr, 10); - fputs("\np = ", stderr); mp_writefile(p.m, stderr, 10); - fputs("\nr = ", stderr); mp_writefile(r, stderr, 10); + fputs("\np = ", stderr); mp_writefile(p, stderr, 10); + fputs("\nq = ", stderr); mp_writefile(q, stderr, 10); fputc('\n', stderr); ok = 0; } mp_drop(m); - mp_drop(r); - pgen_destroy(&p); + mp_drop(q); + mp_drop(p); assert(mparena_count(MPARENA_GLOBAL) == 0); return (ok); } static test_chunk tests[] = { { "pgen", verify, { &type_mp, &type_mp, 0 } }, + { "primep", t_primep, { &type_mp, &type_int, 0 } }, { 0, 0, { 0 } } }; @@ -342,7 +417,6 @@ int main(int argc, char *argv[]) test_run(argc, argv, tests, SRCDIR "/tests/pgen"); return (0); } - #endif /*----- That's all, folks -------------------------------------------------*/