X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/0a5d0bc7225e8abcab0519447a82d87c6d05884a..2685767a6125c1620719c7de6234aedf41857b7e:/mpbarrett.c diff --git a/mpbarrett.c b/mpbarrett.c index 711b7ab..889b268 100644 --- a/mpbarrett.c +++ b/mpbarrett.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mpbarrett.c,v 1.3 1999/12/12 15:08:52 mdw Exp $ + * $Id: mpbarrett.c,v 1.7 2001/04/19 18:25:26 mdw Exp $ * * Barrett modular reduction * @@ -30,6 +30,21 @@ /*----- 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 + * (e.g., replace MP_MODIFY by MP_DEST). + * * Revision 1.3 1999/12/12 15:08:52 mdw * Don't bother shifting %$q$% in @mpbarrett_reduce@, just skip the least * significant digits. @@ -73,7 +88,7 @@ void mpbarrett_create(mpbarrett *mb, mp *m) mp_shrink(m); mb->k = MP_LEN(m); mb->m = MP_COPY(m); - b = mp_create(2 * mb->k + 1); + b = mp_new(2 * mb->k + 1, 0); MPX_ZERO(b->v, b->vl - 1); b->vl[-1] = 1; mp_div(&b, 0, b, m); @@ -105,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) @@ -118,7 +132,8 @@ mp *mpbarrett_reduce(mpbarrett *mb, mp *d, mp *m) if (MP_LEN(m) < k) { m = MP_COPY(m); - MP_DROP(d); + if (d) + MP_DROP(d); return (m); } @@ -130,7 +145,8 @@ mp *mpbarrett_reduce(mpbarrett *mb, mp *d, mp *m) q = mp_mul(MP_NEW, &qq, mb->mu); if (MP_LEN(q) <= k) { m = MP_COPY(m); - MP_DROP(d); + if (d) + MP_DROP(d); return (m); } } @@ -146,12 +162,11 @@ mp *mpbarrett_reduce(mpbarrett *mb, mp *d, mp *m) mvl = m->vl; else mvl = m->v + k + 1; - r = mp_create(k + 1); + r = mp_new(k + 1, (q->f | mb->m->f) & MP_BURN); mpx_umul(r->v, r->vl, q->v + k + 1, q->vl, mb->m->v, mb->m->vl); - r->f = (q->f | mb->m->f) & MP_BURN; - MP_MODIFY(d, k + 1); + 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); @@ -163,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); } @@ -177,35 +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 (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); @@ -214,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 @@ -230,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); @@ -262,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);