X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/ba6e6b64033b1f9de49feccb5c9cd438354481f7..0f00dc4c8eb47e67bc0f148c2dd109f73a451e0a:/math/pgen.c diff --git a/math/pgen.c b/math/pgen.c new file mode 100644 index 0000000..9a822f5 --- /dev/null +++ b/math/pgen.c @@ -0,0 +1,420 @@ +/* -*-c-*- + * + * Prime generation glue + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * 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. + */ + +/*----- Header files ------------------------------------------------------*/ + +#include +#include +#include +#include + +#include "fibrand.h" +#include "grand.h" +#include "mp.h" +#include "mprand.h" +#include "pgen.h" +#include "pfilt.h" +#include "rabin.h" + +/*----- Standard prime filter ---------------------------------------------*/ + +/* --- @pgen_filter@ --- */ + +int pgen_filter(int rq, pgen_event *ev, void *p) +{ + pgen_filterctx *f = p; + int rc = PGEN_FAIL; + + 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; + } + + 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_jump@ --- * + * + * Similar to the standard @pgen_filter@, but jumps in large steps rather + * than small ones. + */ + +int pgen_jump(int rq, pgen_event *ev, void *p) +{ + 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); + } + + while (rc == PGEN_FAIL) + rc = pfilt_jump(&f->f, f->j); + ev->m = MP_COPY(f->f.m); + return (rc); +} + +/*----- Standard prime test -----------------------------------------------*/ + +/* --- @pgen_test@ --- */ + +int pgen_test(int rq, pgen_event *ev, void *p) +{ + rabin *r = p; + int rc = PGEN_ABORT; + + 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; + } + + return (rc); +} + +/*----- The main driver ---------------------------------------------------*/ + +/* --- @pgen@ --- * + * + * 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 + * + * Returns: Pointer to final result, or null. + * + * Use: A generalized prime-number search skeleton. Yes, that's a + * scary number of arguments. + */ + +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) +{ + 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); + } + + /* --- 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; + } + + /* --- If decrementing counters is requested, do that --- */ + + if ((act & A_STEP) && steps) { + ev.steps++; + if (ev.steps == steps) { + act |= A_EVENT | A_ENDSTEP | A_DONE; + rc = PGEN_ABORT; + } + ev.tests = 0; + } + + if ((act & A_TEST) && tests) { + ev.tests++; + if (ev.tests == tests) { + act |= A_ENDTEST | A_ENDSTEP | A_DONE; + rc = PGEN_DONE; + } + } + + /* --- 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; + } + } + + /* --- Close down tester and stepper functions --- */ + + if (act & A_ENDTEST) + test(PGEN_DONE, &ev, tctx); + if (act & A_ENDSTEP) + step(PGEN_DONE, &ev, sctx); + + /* --- Stop the entire test if necessary --- */ + + if (act & A_DONE) + break; + } + + /* --- Tidy up and return --- */ + + if (rc == PGEN_ABORT) { + mp_drop(ev.m); + ev.m = 0; + } + ev.r->ops->destroy(ev.r); + mp_drop(d); + + return (ev.m); +} + +/* --- @pgen_primep@ --- * + * + * Arguments: @mp *p@ = a number to check + * @grand *gr@ = a random number source + * + * Returns: Nonzero if @p@ is really prime. + */ + +int pgen_primep(mp *p, grand *gr) +{ + int i; + rabin r; + mp *x = MP_NEW; + + 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 rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +#include + +static int t_primep(dstr *v) +{ + mp *m = *(mp **)v[0].buf; + int e = *(int *)v[1].buf; + int r; + grand *rng; + int ok = 1; + + 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; + } + + 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, stderr, 10); + fputs("\nq = ", stderr); mp_writefile(q, stderr, 10); + fputc('\n', stderr); + ok = 0; + } + + mp_drop(m); + 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 } } +}; + +int main(int argc, char *argv[]) +{ + sub_init(); + test_run(argc, argv, tests, SRCDIR "/t/pgen"); + return (0); +} +#endif + +/*----- That's all, folks -------------------------------------------------*/