X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/d3409d5ecf2492cff862616de72a580d1a8e8dc0..24e1e862d1cf06eb08af7c44468c1d12fcffc747:/mpmont.c diff --git a/mpmont.c b/mpmont.c index d924009..35bffc2 100644 --- a/mpmont.c +++ b/mpmont.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mpmont.c,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * $Id: mpmont.c,v 1.4 1999/11/21 12:27:06 mdw Exp $ * * Montgomery reduction * @@ -30,6 +30,19 @@ /*----- Revision history --------------------------------------------------* * * $Log: mpmont.c,v $ + * Revision 1.4 1999/11/21 12:27:06 mdw + * Remove a division from the Montgomery setup by calculating + * %$R^2 \bmod m$% first and then %$R \bmod m$% by Montgomery reduction of + * %$R^2$%. + * + * Revision 1.3 1999/11/21 11:35:10 mdw + * Performance improvement: use @mp_sqr@ and @mpmont_reduce@ instead of + * @mpmont_mul@ for squaring in exponentiation. + * + * Revision 1.2 1999/11/19 13:17:26 mdw + * Add extra interface to exponentiation which returns a Montgomerized + * result. + * * Revision 1.1 1999/11/17 18:02:16 mdw * New multiprecision integer arithmetic suite. * @@ -84,15 +97,14 @@ void mpmont_create(mpmont *mm, mp *m) { size_t l = MP_LEN(m); - mp *r = mp_create(l + 1); + mp *r = mp_create(2 * l + 1); mm->shift = l * MPW_BITS; MPX_ZERO(r->v, r->vl - 1); r->vl[-1] = 1; - mm->r = mm->r2 = MP_NEW; - mp_div(0, &mm->r, r, m); - r = mp_sqr(r, mm->r); + mm->r2 = MP_NEW; mp_div(0, &mm->r2, r, m); + mm->r = mpmont_reduce(mm, MP_NEW, mm->r2); MP_DROP(r); } } @@ -134,9 +146,9 @@ mp *mpmont_reduce(mpmont *mm, mp *d, const mp *a) n = MP_LEN(mm->m); if (d == a) - MP_MODIFY(d, 2 * n); + MP_MODIFY(d, 2 * n + 1); else { - MP_MODIFY(d, 2 * n); + MP_MODIFY(d, 2 * n + 1); memcpy(d->v, a->v, MPWS(MP_LEN(a))); memset(d->v + MP_LEN(a), 0, MPWS(MP_LEN(d) - MP_LEN(a))); } @@ -158,6 +170,8 @@ mp *mpmont_reduce(mpmont *mm, mp *d, const mp *a) d->vl -= dv - d->v; MP_SHRINK(d); d->f = a->f & MP_BURN; + if (MP_CMP(d, >=, mm->m)) + d = mp_sub(d, d, mm->m); return (d); } @@ -221,43 +235,68 @@ mp *mpmont_mul(mpmont *mm, mp *d, const mp *a, const mp *b) d->vl -= dv - d->v; MP_SHRINK(d); d->f = (a->f | b->f) & MP_BURN; + if (MP_CMP(d, >=, mm->m)) + d = mp_sub(d, d, mm->m); return (d); } -/* --- @mpmont_exp@ --- * +/* --- @mpmont_expr@ --- * * * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context * @const mp *a@ = base * @const mp *e@ = exponent * - * Returns: Result, %$a^e \bmod m$%. + * Returns: Result, %$a^e R \bmod m$%. */ -mp *mpmont_exp(mpmont *mm, const mp *a, const mp *e) +mp *mpmont_expr(mpmont *mm, const mp *a, const mp *e) { mpscan sc; mp *ar = mpmont_mul(mm, MP_NEW, a, mm->r2); mp *d = MP_COPY(mm->r); + mp *spare = MP_NEW; mp_scan(&sc, e); if (MP_STEP(&sc)) { + size_t sq = 0; for (;;) { mp *dd; if (MP_BIT(&sc)) { - dd = mpmont_mul(mm, MP_NEW, d, ar); - MP_DROP(d); - d = dd; + while (sq) { + dd = mp_sqr(spare, ar); + dd = mpmont_reduce(mm, dd, dd); + spare = ar; ar = dd; + sq--; + } + dd = mpmont_mul(mm, spare, d, ar); + spare = d; d = dd; } + sq++; if (!MP_STEP(&sc)) break; - dd = mpmont_mul(mm, MP_NEW, ar, ar); - MP_DROP(ar); - ar = dd; } } MP_DROP(ar); - return (mpmont_reduce(mm, d, d)); + if (spare != MP_NEW) + MP_DROP(spare); + return (d); +} + +/* --- @mpmont_exp@ --- * + * + * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context + * @const mp *a@ = base + * @const mp *e@ = exponent + * + * Returns: Result, %$a^e \bmod m$%. + */ + +mp *mpmont_exp(mpmont *mm, const mp *a, const mp *e) +{ + mp *d = mpmont_expr(mm, a, e); + d = mpmont_reduce(mm, d, d); + return (d); } /*----- Test rig ----------------------------------------------------------*/ @@ -288,7 +327,7 @@ static int tcreate(dstr *v) fputs("\n*** bad r", stderr); fputs("\nm = ", stderr); mp_writefile(m, stderr, 10); fputs("\nexpected ", stderr); mp_writefile(r, stderr, 10); - fputs("\n found ", stderr); mp_writefile(r, stderr, 10); + fputs("\n found ", stderr); mp_writefile(mm.r, stderr, 10); fputc('\n', stderr); ok = 0; } @@ -297,7 +336,7 @@ static int tcreate(dstr *v) fputs("\n*** bad r2", stderr); fputs("\nm = ", stderr); mp_writefile(m, stderr, 10); fputs("\nexpected ", stderr); mp_writefile(r2, stderr, 10); - fputs("\n found ", stderr); mp_writefile(r2, stderr, 10); + fputs("\n found ", stderr); mp_writefile(mm.r2, stderr, 10); fputc('\n', stderr); ok = 0; } @@ -316,51 +355,53 @@ static int tmul(dstr *v) mp *a = *(mp **)v[1].buf; mp *b = *(mp **)v[2].buf; mp *r = *(mp **)v[3].buf; - mp *mr, *qr; int ok = 1; mpmont mm; mpmont_create(&mm, m); { + mp *qr = mp_mul(MP_NEW, a, b); + mp_div(0, &qr, qr, m); + + if (MP_CMP(qr, !=, r)) { + fputs("\n*** classical modmul failed", stderr); + fputs("\n m = ", stderr); mp_writefile(m, stderr, 10); + fputs("\n a = ", stderr); mp_writefile(a, stderr, 10); + fputs("\n b = ", stderr); mp_writefile(b, stderr, 10); + fputs("\n r = ", stderr); mp_writefile(r, stderr, 10); + fputs("\nqr = ", stderr); mp_writefile(qr, stderr, 10); + fputc('\n', stderr); + ok = 0; + } + + mp_drop(qr); + } + + { mp *ar = mpmont_mul(&mm, MP_NEW, a, mm.r2); mp *br = mpmont_mul(&mm, MP_NEW, b, mm.r2); - mr = mpmont_mul(&mm, MP_NEW, ar, br); + mp *mr = mpmont_mul(&mm, MP_NEW, ar, br); mr = mpmont_reduce(&mm, mr, mr); + if (MP_CMP(mr, !=, r)) { + fputs("\n*** montgomery modmul failed", stderr); + fputs("\n m = ", stderr); mp_writefile(m, stderr, 10); + fputs("\n a = ", stderr); mp_writefile(a, stderr, 10); + fputs("\n b = ", 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(ar); MP_DROP(br); + mp_drop(mr); } - qr = mp_mul(MP_NEW, a, b); - mp_div(0, &qr, qr, m); - - if (MP_CMP(qr, !=, r)) { - fputs("\n*** classical modmul failed", stderr); - fputs("\n m = ", stderr); mp_writefile(m, stderr, 10); - fputs("\n a = ", stderr); mp_writefile(a, stderr, 10); - fputs("\n b = ", stderr); mp_writefile(b, stderr, 10); - fputs("\n r = ", stderr); mp_writefile(r, stderr, 10); - fputs("\nqr = ", stderr); mp_writefile(qr, stderr, 10); - fputc('\n', stderr); - ok = 0; - } - - if (MP_CMP(mr, !=, r)) { - fputs("\n*** montgomery modmul failed", stderr); - fputs("\n m = ", stderr); mp_writefile(m, stderr, 10); - fputs("\n a = ", stderr); mp_writefile(a, stderr, 10); - fputs("\n b = ", 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); - MP_DROP(qr); mpmont_destroy(&mm); return ok; }