X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/d34decd2b2b88240cf4ca68a2a5feb7bf36de6e7..2685767a6125c1620719c7de6234aedf41857b7e:/mpbarrett.c diff --git a/mpbarrett.c b/mpbarrett.c index aaa1bb7..889b268 100644 --- a/mpbarrett.c +++ b/mpbarrett.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mpbarrett.c,v 1.4 2000/06/17 11:45:09 mdw Exp $ + * $Id: mpbarrett.c,v 1.7 2001/04/19 18:25:26 mdw Exp $ * * Barrett modular reduction * @@ -30,6 +30,16 @@ /*----- Revision history --------------------------------------------------* * * $Log: mpbarrett.c,v $ + * Revision 1.7 2001/04/19 18:25:26 mdw + * Use sliding-window exponentiation. + * + * Revision 1.6 2000/10/08 12:03:44 mdw + * (mpbarrett_reduce): Cope with negative numbers. + * + * Revision 1.5 2000/07/29 17:04:33 mdw + * Change to use left-to-right bitwise exponentiation. This will improve + * performance when the base is small. + * * Revision 1.4 2000/06/17 11:45:09 mdw * Major memory management overhaul. Added arena support. Use the secure * arena for secret integers. Replace and improve the MP management macros @@ -110,8 +120,7 @@ void mpbarrett_destroy(mpbarrett *mb) * Returns: The residue of @m@ modulo the number in the reduction * context. * - * Use: Performs an efficient modular reduction. The argument is - * assumed to be positive. + * Use: Performs an efficient modular reduction. */ mp *mpbarrett_reduce(mpbarrett *mb, mp *d, mp *m) @@ -157,7 +166,7 @@ mp *mpbarrett_reduce(mpbarrett *mb, mp *d, mp *m) mpx_umul(r->v, r->vl, q->v + k + 1, q->vl, mb->m->v, mb->m->vl); MP_DEST(d, k + 1, r->f); mpx_usub(d->v, d->vl, m->v, mvl, r->v, r->vl); - d->f = (m->f | r->f) & MP_BURN; + d->f = (m->f | r->f) & (MP_BURN | MP_NEG); MP_DROP(r); MP_DROP(q); MP_DROP(m); @@ -169,6 +178,13 @@ mp *mpbarrett_reduce(mpbarrett *mb, mp *d, mp *m) while (MPX_UCMP(d->v, d->vl, >=, mb->m->v, mb->m->vl)) mpx_usub(d->v, d->vl, d->v, d->vl, mb->m->v, mb->m->vl); + /* --- Fix up the sign --- */ + + if (d->f & MP_NEG) { + mpx_usub(d->v, d->vl, mb->m->v, mb->m->vl, d->v, d->vl); + d->f &= ~MP_NEG; + } + MP_SHRINK(d); return (d); } @@ -183,37 +199,62 @@ mp *mpbarrett_reduce(mpbarrett *mb, mp *d, mp *m) * Returns: Result, %$a^e \bmod m$%. */ -mp *mpbarrett_exp(mpbarrett *mb, mp *d, mp *a, mp *e) +#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 = MP_NEW; + mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW; + unsigned sq = 0; a = MP_COPY(a); - mp_scan(&sc, e); - if (e->f & MP_BURN) - spare = MP_NEWSEC; - if (MP_STEP(&sc)) { - size_t sq = 0; + 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 (;;) { - mp *dd; - if (MP_BIT(&sc)) { - while (sq) { - dd = mp_sqr(spare, a); - dd = mpbarrett_reduce(mb, dd, dd); - spare = a; a = dd; - sq--; - } - dd = mp_mul(spare, x, a); - dd = mpbarrett_reduce(mb, dd, dd); - spare = x; x = dd; - } + if (!MP_RSTEP(&sc)) + goto done; + if (MP_RBIT(&sc)) + break; sq++; - if (!MP_STEP(&sc)) - break; } } + /* --- 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); @@ -222,6 +263,130 @@ mp *mpbarrett_exp(mpbarrett *mb, mp *d, mp *a, mp *e) 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 @@ -238,7 +403,7 @@ static int vmod(dstr *v) mpbarrett_create(&mb, n); s = mpbarrett_reduce(&mb, MP_NEW, x); - if (MP_CMP(s, !=, r)) { + if (!MP_EQ(s, r)) { fputs("\n*** barrett reduction failure\n", stderr); fputs("x = ", stderr); mp_writefile(x, stderr, 10); fputc('\n', stderr); fputs("n = ", stderr); mp_writefile(n, stderr, 10); fputc('\n', stderr); @@ -270,7 +435,7 @@ static int vexp(dstr *v) mr = mpbarrett_exp(&mb, MP_NEW, a, b); - if (MP_CMP(mr, !=, r)) { + 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);