X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/1ba69a1e60cc3e3b8e1ddff0689ae537e306299b..34e4f738bcba58e6d8c4cabbb0b3232a65b42a9d:/mpbarrett.c diff --git a/mpbarrett.c b/mpbarrett.c index 889b268..934097d 100644 --- a/mpbarrett.c +++ b/mpbarrett.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mpbarrett.c,v 1.7 2001/04/19 18:25:26 mdw Exp $ + * $Id: mpbarrett.c,v 1.9 2004/04/01 12:50:09 mdw Exp $ * * Barrett modular reduction * @@ -30,6 +30,16 @@ /*----- Revision history --------------------------------------------------* * * $Log: mpbarrett.c,v $ + * Revision 1.9 2004/04/01 12:50:09 mdw + * Add cyclic group abstraction, with test code. Separate off exponentation + * functions for better static linking. Fix a buttload of bugs on the way. + * Generally ensure that negative exponents do inversion correctly. Add + * table of standard prime-field subgroups. (Binary field subgroups are + * currently unimplemented but easy to add if anyone ever finds a good one.) + * + * Revision 1.8 2001/06/16 13:00:20 mdw + * Use the generic exponentiation functions. + * * Revision 1.7 2001/04/19 18:25:26 mdw * Use sliding-window exponentiation. * @@ -189,204 +199,6 @@ mp *mpbarrett_reduce(mpbarrett *mb, mp *d, mp *m) return (d); } -/* --- @mpbarrett_exp@ --- * - * - * Arguments: @mpbarrett *mb@ = pointer to Barrett reduction context - * @mp *d@ = fake destination - * @mp *a@ = base - * @mp *e@ = exponent - * - * Returns: Result, %$a^e \bmod m$%. - */ - -#define WINSZ 5 -#define TABSZ (1 << (WINSZ - 1)) - -#define THRESH (((MPW_BITS / WINSZ) << 2) + 1) - -static mp *exp_simple(mpbarrett *mb, mp *d, mp *a, mp *e) -{ - mpscan sc; - mp *x = MP_ONE; - mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW; - unsigned sq = 0; - - a = MP_COPY(a); - mp_rscan(&sc, e); - if (!MP_RSTEP(&sc)) - goto exit; - while (!MP_RBIT(&sc)) - MP_RSTEP(&sc); - - /* --- Do the main body of the work --- */ - - for (;;) { - sq++; - while (sq) { - mp *y; - y = mp_sqr(spare, x); - y = mpbarrett_reduce(mb, y, y); - spare = x; x = y; - sq--; - } - { - mp *y = mp_mul(spare, x, a); - y = mpbarrett_reduce(mb, y, y); - spare = x; x = y; - } - for (;;) { - if (!MP_RSTEP(&sc)) - goto done; - if (MP_RBIT(&sc)) - break; - sq++; - } - } - - /* --- Do a final round of squaring --- */ - -done: - while (sq) { - mp *y; - y = mp_sqr(spare, x); - y = mpbarrett_reduce(mb, y, y); - spare = x; x = y; - sq--; - } - -exit: - MP_DROP(a); - if (spare != MP_NEW) - MP_DROP(spare); - if (d != MP_NEW) - MP_DROP(d); - return (x); -} - -mp *mpbarrett_exp(mpbarrett *mb, mp *d, mp *a, mp *e) -{ - mp **tab; - mp *a2; - mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW; - mp *x = MP_ONE; - unsigned i, sq = 0; - mpscan sc; - - /* --- Do we bother? --- */ - - MP_SHRINK(e); - if (MP_LEN(e) == 0) - goto exit; - if (MP_LEN(e) < THRESH) { - x->ref--; - return (exp_simple(mb, d, a, e)); - } - - /* --- Do the precomputation --- */ - - a2 = mp_sqr(MP_NEW, a); - a2 = mpbarrett_reduce(mb, a2, a2); - tab = xmalloc(TABSZ * sizeof(mp *)); - tab[0] = MP_COPY(a); - for (i = 1; i < TABSZ; i++) { - mp *x = mp_mul(MP_NEW, tab[i - 1], a2); - tab[i] = mpbarrett_reduce(mb, x, x); - } - mp_drop(a2); - mp_rscan(&sc, e); - - /* --- Skip top-end zero bits --- * - * - * If the initial step worked, there must be a set bit somewhere, so keep - * stepping until I find it. - */ - - MP_RSTEP(&sc); - while (!MP_RBIT(&sc)) - MP_RSTEP(&sc); - - /* --- Now for the main work --- */ - - for (;;) { - unsigned l = 0; - unsigned z = 0; - - /* --- The next bit is set, so read a window index --- * - * - * Reset @i@ to zero and increment @sq@. Then, until either I read - * @WINSZ@ bits or I run out of bits, scan in a bit: if it's clear, bump - * the @z@ counter; if it's set, push a set bit into @i@, shift it over - * by @z@ bits, bump @sq@ by @z + 1@ and clear @z@. By the end of this - * palaver, @i@ is an index to the precomputed value in @tab@. - */ - - i = 0; - sq++; - for (;;) { - l++; - if (l >= WINSZ || !MP_RSTEP(&sc)) - break; - if (!MP_RBIT(&sc)) - z++; - else { - i = ((i << 1) | 1) << z; - sq += z + 1; - z = 0; - } - } - - /* --- Do the squaring --- * - * - * Remember that @sq@ carries over from the zero-skipping stuff below. - */ - - while (sq) { - mp *y; - y = mp_sqr(spare, x); - y = mpbarrett_reduce(mb, y, y); - spare = x; x = y; - sq--; - } - - /* --- Do the multiply --- */ - - { mp *y = mp_mul(spare, x, tab[i]); spare = x; - x = mpbarrett_reduce(mb, y, y); } - - /* --- Now grind along through the rest of the bits --- */ - - sq = z; - for (;;) { - if (!MP_RSTEP(&sc)) - goto done; - if (MP_RBIT(&sc)) - break; - sq++; - } - } - - /* --- Do a final round of squaring --- */ - -done: - while (sq) { - mp *y; - y = mp_sqr(spare, x); - y = mpbarrett_reduce(mb, y, y); - spare = x; x = y; - sq--; - } - - /* --- Done --- */ - - for (i = 0; i < TABSZ; i++) - mp_drop(tab[i]); - xfree(tab); -exit: - mp_drop(d); - mp_drop(spare); - return (x); -} - /*----- Test rig ----------------------------------------------------------*/ #ifdef TEST_RIG @@ -421,44 +233,8 @@ static int vmod(dstr *v) return (ok); } -static int vexp(dstr *v) -{ - mp *m = *(mp **)v[0].buf; - mp *a = *(mp **)v[1].buf; - mp *b = *(mp **)v[2].buf; - mp *r = *(mp **)v[3].buf; - mp *mr; - int ok = 1; - - mpbarrett mb; - mpbarrett_create(&mb, m); - - mr = mpbarrett_exp(&mb, MP_NEW, a, b); - - if (!MP_EQ(mr, r)) { - fputs("\n*** barrett modexp failed", stderr); - fputs("\n m = ", stderr); mp_writefile(m, stderr, 10); - fputs("\n a = ", stderr); mp_writefile(a, stderr, 10); - fputs("\n e = ", stderr); mp_writefile(b, stderr, 10); - fputs("\n r = ", stderr); mp_writefile(r, stderr, 10); - fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10); - fputc('\n', stderr); - ok = 0; - } - - mp_drop(m); - mp_drop(a); - mp_drop(b); - mp_drop(r); - mp_drop(mr); - mpbarrett_destroy(&mb); - assert(mparena_count(MPARENA_GLOBAL) == 0); - return ok; -} - static test_chunk tests[] = { { "mpbarrett-reduce", vmod, { &type_mp, &type_mp, &type_mp, 0 } }, - { "mpbarrett-exp", vexp, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } }, { 0, 0, { 0 } } };