X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/f5f35081384158156cd62ef93e51c1afce1b9aa9..c9d4c30b618d250d16676a9db4ef225cee9ec77c:/mpmont.c diff --git a/mpmont.c b/mpmont.c index a5b5e66..be7d00e 100644 --- a/mpmont.c +++ b/mpmont.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mpmont.c,v 1.7 1999/12/11 01:51:14 mdw Exp $ + * $Id: mpmont.c,v 1.10 2000/07/29 17:05:43 mdw Exp $ * * Montgomery reduction * @@ -30,6 +30,19 @@ /*----- Revision history --------------------------------------------------* * * $Log: mpmont.c,v $ + * Revision 1.10 2000/07/29 17:05:43 mdw + * (mpmont_expr): Use sliding window exponentiation, with a drop-through + * for small exponents to use a simple left-to-right bitwise routine. This + * can reduce modexp times by up to a quarter. + * + * Revision 1.9 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.8 1999/12/22 15:55:00 mdw + * Adjust Karatsuba parameters. + * * Revision 1.7 1999/12/11 01:51:14 mdw * Use a Karatsuba-based reduction for large moduli. * @@ -102,7 +115,7 @@ void mpmont_create(mpmont *mm, mp *m) void mpmont_create(mpmont *mm, mp *m) { size_t n = MP_LEN(m); - mp *r2 = mp_create(2 * n + 1); + mp *r2 = mp_new(2 * n + 1, 0); mp r; /* --- Validate the arguments --- */ @@ -182,7 +195,7 @@ mp *mpmont_reduce(mpmont *mm, mp *d, mp *a) /* --- Check for serious Karatsuba reduction --- */ - if (n > KARATSUBA_CUTOFF * 2) { + if (n > KARATSUBA_CUTOFF * 3) { mp al; mpw *vl; mp *u; @@ -210,13 +223,12 @@ mp *mpmont_reduce(mpmont *mm, mp *d, mp *a) /* --- Initial conditioning of the arguments --- */ - if (d == a) - MP_MODIFY(d, 2 * n + 1); - else { - MP_MODIFY(d, 2 * n + 1); - MPX_COPY(d->v, d->vl, a->v, a->vl); - } - + a = MP_COPY(a); + if (d) + MP_DROP(d); + d = a; + MP_DEST(d, 2 * n + 1, a->f); + dv = d->v; dvl = d->vl; mv = mm->m->v; mvl = mm->m->vl; @@ -232,7 +244,6 @@ mp *mpmont_reduce(mpmont *mm, mp *d, mp *a) /* --- Wrap everything up --- */ - d->f = a->f & MP_BURN; memmove(d->v, d->v + n, MPWS(MP_LEN(d) - n)); d->vl -= n; if (MP_CMP(d, >=, mm->m)) @@ -265,7 +276,7 @@ mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b) mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b) { - if (mm->n > KARATSUBA_CUTOFF * 2) { + if (mm->n > KARATSUBA_CUTOFF * 3) { d = mp_mul(d, a, b); d = mpmont_reduce(mm, d, d); } else { @@ -286,7 +297,7 @@ mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b) a = MP_COPY(a); b = MP_COPY(b); - MP_MODIFY(d, 2 * n + 1); + MP_DEST(d, 2 * n + 1, a->f | b->f | MP_UNDEF); dv = d->v; dvl = d->vl; MPX_ZERO(dv, dvl); av = a->v; avl = a->vl; @@ -343,35 +354,63 @@ mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b) * Returns: Result, %$a^e R \bmod m$%. */ -mp *mpmont_expr(mpmont *mm, 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(mpmont *mm, mp *d, mp *a, mp *e) { mpscan sc; - mp *ar = mpmont_mul(mm, MP_NEW, a, mm->r2); + mp *ar; mp *x = MP_COPY(mm->r); - mp *spare = MP_NEW; - - mp_scan(&sc, e); - - if (MP_STEP(&sc)) { - size_t sq = 0; + mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW; + unsigned 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 --- */ + + ar = mpmont_mul(mm, MP_NEW, a, mm->r2); + for (;;) { + sq++; + while (sq) { + mp *y; + y = mp_sqr(spare, x); + y = mpmont_reduce(mm, y, y); + spare = x; x = y; + sq--; + } + { mp *y = mpmont_mul(mm, spare, x, ar); spare = x; x = y; } + sq = 0; for (;;) { - mp *dd; - if (MP_BIT(&sc)) { - while (sq) { - dd = mp_sqr(spare, ar); - dd = mpmont_reduce(mm, dd, dd); - spare = ar; ar = dd; - sq--; - } - dd = mpmont_mul(mm, spare, x, ar); - spare = x; x = dd; - } - sq++; - if (!MP_STEP(&sc)) + 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 = mpmont_reduce(mm, y, y); + spare = x; x = y; + sq--; + } + + /* --- Done --- */ + MP_DROP(ar); +exit: if (spare != MP_NEW) MP_DROP(spare); if (d != MP_NEW) @@ -379,6 +418,131 @@ mp *mpmont_expr(mpmont *mm, mp *d, mp *a, mp *e) return (x); } +mp *mpmont_expr(mpmont *mm, mp *d, mp *a, mp *e) +{ + mp **tab; + mp *ar, *a2; + mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW; + mp *x = MP_COPY(mm->r); + 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(mm, d, a, e)); + } + + /* --- Do the precomputation --- */ + + ar = mpmont_mul(mm, MP_NEW, a, mm->r2); + a2 = mp_sqr(MP_NEW, ar); + a2 = mpmont_reduce(mm, a2, a2); + tab = xmalloc(TABSZ * sizeof(mp *)); + tab[0] = ar; + for (i = 1; i < TABSZ; i++) + tab[i] = mpmont_mul(mm, MP_NEW, tab[i - 1], a2); + 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 = mpmont_reduce(mm, y, y); + spare = x; x = y; + sq--; + } + + /* --- Do the multiply --- */ + + { mp *y = mpmont_mul(mm, spare, x, tab[i]); spare = x; x = 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 = mpmont_reduce(mm, y, y); + spare = x; x = y; + sq--; + } + + /* --- Done --- */ + + for (i = 0; i < TABSZ; i++) + mp_drop(tab[i]); + xfree(tab); +exit: + if (d != MP_NEW) + mp_drop(d); + if (spare) + mp_drop(spare); + return (x); +} + /* --- @mpmont_exp@ --- * * * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context