X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/17ad212e749b467593355532eeb475374277463e..ef5f48103e83977bda6ef4d7d1aacbb66a629b10:/mpmont.c diff --git a/mpmont.c b/mpmont.c index 4885151..d061934 100644 --- a/mpmont.c +++ b/mpmont.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mpmont.c,v 1.2 1999/11/19 13:17:26 mdw Exp $ + * $Id: mpmont.c,v 1.6 1999/12/10 23:18:39 mdw Exp $ * * Montgomery reduction * @@ -30,6 +30,22 @@ /*----- Revision history --------------------------------------------------* * * $Log: mpmont.c,v $ + * Revision 1.6 1999/12/10 23:18:39 mdw + * Change interface for suggested destinations. + * + * Revision 1.5 1999/11/22 13:58:40 mdw + * Add an option to disable Montgomery reduction, so that performance + * comparisons can be done. + * + * 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. @@ -44,6 +60,16 @@ #include "mp.h" #include "mpmont.h" +/*----- Tweakables --------------------------------------------------------*/ + +/* --- @MPMONT_DISABLE@ --- * + * + * Replace all the clever Montgomery reduction with good old-fashioned long + * division. + */ + +/* #define MPMONT_DISABLE */ + /*----- Main code ---------------------------------------------------------*/ /* --- @mpmont_create@ --- * @@ -54,10 +80,29 @@ * Returns: --- * * Use: Initializes a Montgomery reduction context ready for use. + * The argument @m@ must be a positive odd integer. */ +#ifdef MPMONT_DISABLE + void mpmont_create(mpmont *mm, mp *m) { + mp_shrink(m); + mm->m = MP_COPY(m); + mm->r = MP_ONE; + mm->r2 = MP_ONE; +} + +#else + +void mpmont_create(mpmont *mm, mp *m) +{ + /* --- Validate the arguments --- */ + + assert(((void)"Montgomery modulus must be positive", + (m->f & MP_NEG) == 0)); + assert(((void)"Montgomery modulus must be odd", m->v[0] & 1)); + /* --- Take a copy of the modulus --- */ mp_shrink(m); @@ -71,7 +116,7 @@ void mpmont_create(mpmont *mm, mp *m) { mpw av[2] = { 0, 1 }; mp a, b; - mp *i; + mp *i = MP_NEW; mpw mi; mp_build(&a, av, av + 2); @@ -88,19 +133,20 @@ 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); } } +#endif + /* --- @mpmont_destroy@ --- * * * Arguments: @mpmont *mm@ = pointer to a Montgomery reduction context @@ -122,15 +168,25 @@ void mpmont_destroy(mpmont *mm) * * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context * @mp *d@ = destination - * @const mp *a@ = source, assumed positive + * @mp *a@ = source, assumed positive * * Returns: Result, %$a R^{-1} \bmod m$%. */ -mp *mpmont_reduce(mpmont *mm, mp *d, const mp *a) +#ifdef MPMONT_DISABLE + +mp *mpmont_reduce(mpmont *mm, mp *d, mp *a) +{ + mp_div(0, &d, a, mm->m); + return (d); +} + +#else + +mp *mpmont_reduce(mpmont *mm, mp *d, mp *a) { mpw *dv, *dvl; - const mpw *mv, *mvl; + mpw *mv, *mvl; size_t n; /* --- Initial conditioning of the arguments --- */ @@ -138,9 +194,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))); } @@ -167,85 +223,111 @@ mp *mpmont_reduce(mpmont *mm, mp *d, const mp *a) return (d); } +#endif + /* --- @mpmont_mul@ --- * * * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context * @mp *d@ = destination - * @const mp *a, *b@ = sources, assumed positive + * @mp *a, *b@ = sources, assumed positive * * Returns: Result, %$a b R^{-1} \bmod m$%. */ -mp *mpmont_mul(mpmont *mm, mp *d, const mp *a, const mp *b) +#ifdef MPMONT_DISABLE + +mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b) { - mpw *dv, *dvl; - const mpw *av, *avl; - const mpw *bv, *bvl; - const mpw *mv, *mvl; - mpw y; - size_t n, i; + d = mp_mul(d, a, b); + mp_div(0, &d, d, mm->m); + return (d); +} - /* --- Initial conditioning of the arguments --- */ +#else - if (MP_LEN(a) > MP_LEN(b)) { - const mp *t = a; a = b; b = t; - } - n = MP_LEN(mm->m); - - MP_MODIFY(d, 2 * n + 1); - dv = d->v; dvl = d->vl; - MPX_ZERO(dv, dvl); - av = a->v; avl = a->vl; - bv = b->v; bvl = b->vl; - mv = mm->m->v; mvl = mm->m->vl; - y = *bv; +mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b) +{ + if (MP_LEN(a) > KARATSUBA_CUTOFF && MP_LEN(b) > KARATSUBA_CUTOFF) { + d = mp_mul(d, a, b); + d = mpmont_reduce(mm, d, d); + } else { + mpw *dv, *dvl; + mpw *av, *avl; + mpw *bv, *bvl; + mpw *mv, *mvl; + mpw y; + size_t n, i; + + /* --- Initial conditioning of the arguments --- */ + + if (MP_LEN(a) > MP_LEN(b)) { + mp *t = a; a = b; b = t; + } + n = MP_LEN(mm->m); + + a = MP_COPY(a); + b = MP_COPY(b); + MP_MODIFY(d, 2 * n + 1); + dv = d->v; dvl = d->vl; + MPX_ZERO(dv, dvl); + av = a->v; avl = a->vl; + bv = b->v; bvl = b->vl; + mv = mm->m->v; mvl = mm->m->vl; + y = *bv; + + /* --- Montgomery multiplication phase --- */ + + i = 0; + while (i < n && av < avl) { + mpw x = *av++; + mpw u = MPW((*dv + x * y) * mm->mi); + MPX_UMLAN(dv, dvl, bv, bvl, x); + MPX_UMLAN(dv, dvl, mv, mvl, u); + dv++; + i++; + } - /* --- Montgomery multiplication phase --- */ + /* --- Simpler Montgomery reduction phase --- */ - i = 0; - while (i < n && av < avl) { - mpw x = *av++; - mpw u = MPW((*dv + x * y) * mm->mi); - MPX_UMLAN(dv, dvl, bv, bvl, x); - MPX_UMLAN(dv, dvl, mv, mvl, u); - dv++; - i++; - } + while (i < n) { + mpw u = MPW(*dv * mm->mi); + MPX_UMLAN(dv, dvl, mv, mvl, u); + dv++; + i++; + } - /* --- Simpler Montgomery reduction phase --- */ + /* --- Done --- */ - while (i < n) { - mpw u = MPW(*dv * mm->mi); - MPX_UMLAN(dv, dvl, mv, mvl, u); - dv++; - i++; + memmove(d->v, dv, MPWS(dvl - dv)); + 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); + MP_DROP(a); + MP_DROP(b); } - /* --- Done --- */ - - memmove(d->v, dv, MPWS(dvl - dv)); - 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); } +#endif + /* --- @mpmont_expr@ --- * * * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context - * @const mp *a@ = base - * @const mp *e@ = exponent + * @mp *d@ = fake destination + * @mp *a@ = base + * @mp *e@ = exponent * * Returns: Result, %$a^e R \bmod m$%. */ -mp *mpmont_expr(mpmont *mm, const mp *a, const mp *e) +mp *mpmont_expr(mpmont *mm, mp *d, mp *a, mp *e) { mpscan sc; mp *ar = mpmont_mul(mm, MP_NEW, a, mm->r2); - mp *d = MP_COPY(mm->r); + mp *x = MP_COPY(mm->r); mp *spare = MP_NEW; mp_scan(&sc, e); @@ -256,12 +338,13 @@ mp *mpmont_expr(mpmont *mm, const mp *a, const mp *e) mp *dd; if (MP_BIT(&sc)) { while (sq) { - dd = mpmont_mul(mm, spare, ar, ar); + 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; + dd = mpmont_mul(mm, spare, x, ar); + spare = x; x = dd; } sq++; if (!MP_STEP(&sc)) @@ -271,21 +354,24 @@ mp *mpmont_expr(mpmont *mm, const mp *a, const mp *e) MP_DROP(ar); if (spare != MP_NEW) MP_DROP(spare); - return (d); + if (d != MP_NEW) + MP_DROP(d); + return (x); } /* --- @mpmont_exp@ --- * * * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context - * @const mp *a@ = base - * @const mp *e@ = exponent + * @mp *d@ = fake destination + * @mp *a@ = base + * @mp *e@ = exponent * * Returns: Result, %$a^e \bmod m$%. */ -mp *mpmont_exp(mpmont *mm, const mp *a, const mp *e) +mp *mpmont_exp(mpmont *mm, mp *d, mp *a, mp *e) { - mp *d = mpmont_expr(mm, a, e); + d = mpmont_expr(mm, d, a, e); d = mpmont_reduce(mm, d, d); return (d); } @@ -337,6 +423,7 @@ static int tcreate(dstr *v) MP_DROP(r); MP_DROP(r2); mpmont_destroy(&mm); + assert(mparena_count(MPARENA_GLOBAL) == 0); return (ok); } @@ -346,52 +433,55 @@ 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); + assert(mparena_count(MPARENA_GLOBAL) == 0); return ok; } @@ -407,7 +497,7 @@ static int texp(dstr *v) mpmont mm; mpmont_create(&mm, m); - mr = mpmont_exp(&mm, a, b); + mr = mpmont_exp(&mm, MP_NEW, a, b); if (MP_CMP(mr, !=, r)) { fputs("\n*** montgomery modexp failed", stderr); @@ -426,14 +516,15 @@ static int texp(dstr *v) MP_DROP(r); MP_DROP(mr); mpmont_destroy(&mm); + assert(mparena_count(MPARENA_GLOBAL) == 0); return ok; } static test_chunk tests[] = { - { "create", tcreate, { &type_mp, &type_mp, &type_mp, &type_mp } }, - { "mul", tmul, { &type_mp, &type_mp, &type_mp, &type_mp } }, - { "exp", texp, { &type_mp, &type_mp, &type_mp, &type_mp } }, + { "create", tcreate, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } }, + { "mul", tmul, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } }, + { "exp", texp, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } }, { 0, 0, { 0 } }, };