X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/17ad212e749b467593355532eeb475374277463e..4640a0dd947f422cfda780c62da342792af27b7b:/mpmont.c diff --git a/mpmont.c b/mpmont.c index 4885151..2371de8 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.15 2001/06/16 13:00:20 mdw Exp $ * * Montgomery reduction * @@ -30,6 +30,54 @@ /*----- Revision history --------------------------------------------------* * * $Log: mpmont.c,v $ + * Revision 1.15 2001/06/16 13:00:20 mdw + * Use the generic exponentiation functions. + * + * Revision 1.14 2001/02/22 09:04:26 mdw + * Cosmetic fix. + * + * Revision 1.13 2001/02/03 12:00:29 mdw + * Now @mp_drop@ checks its argument is non-NULL before attempting to free + * it. Note that the macro version @MP_DROP@ doesn't do this. + * + * Revision 1.12 2000/10/08 15:48:35 mdw + * Rename Karatsuba constants now that we have @gfx_kmul@ too. + * + * Revision 1.11 2000/10/08 12:04:27 mdw + * (mpmont_reduce, mpmont_mul): Cope with negative numbers. + * + * 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. + * + * 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. @@ -43,8 +91,19 @@ #include "mp.h" #include "mpmont.h" +#include "mpmont-exp.h" + +/*----- Tweakables --------------------------------------------------------*/ + +/* --- @MPMONT_DISABLE@ --- * + * + * Replace all the clever Montgomery reduction with good old-fashioned long + * division. + */ + +/* #define MPMONT_DISABLE */ -/*----- Main code ---------------------------------------------------------*/ +/*----- Reduction and multiplication --------------------------------------*/ /* --- @mpmont_create@ --- * * @@ -54,53 +113,62 @@ * 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; + mm->mi = MP_ONE; +} + +#else + void mpmont_create(mpmont *mm, mp *m) { + size_t n = MP_LEN(m); + mp *r2 = mp_new(2 * n + 1, 0); + mp r; + + /* --- 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); mm->m = MP_COPY(m); - /* --- Find the magic value @mi@ --- * - * - * This is a slightly grungy way of solving the problem, but it does work. - */ + /* --- Determine %$R^2$% --- */ - { - mpw av[2] = { 0, 1 }; - mp a, b; - mp *i; - mpw mi; + mm->n = n; + MPX_ZERO(r2->v, r2->vl - 1); + r2->vl[-1] = 1; - mp_build(&a, av, av + 2); - mp_build(&b, m->v, m->v + 1); - mp_gcd(0, 0, &i, &a, &b); - mi = i->v[0]; - if (!(i->f & MP_NEG)) - mi = MPW(-mi); - mm->mi = mi; - MP_DROP(i); - } + /* --- Find the magic value @mi@ --- */ + + mp_build(&r, r2->v + n, r2->vl); + mm->mi = MP_NEW; + mp_gcd(0, 0, &mm->mi, &r, m); + mm->mi = mp_sub(mm->mi, &r, mm->mi); /* --- Discover the values %$R \bmod m$% and %$R^2 \bmod m$% --- */ - { - size_t l = MP_LEN(m); - mp *r = mp_create(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); - mp_div(0, &mm->r2, r, m); - MP_DROP(r); - } + mm->r2 = MP_NEW; + mp_div(0, &mm->r2, r2, m); + mm->r = mpmont_reduce(mm, MP_NEW, mm->r2); + MP_DROP(r2); } +#endif + /* --- @mpmont_destroy@ --- * * * Arguments: @mpmont *mm@ = pointer to a Montgomery reduction context @@ -116,176 +184,232 @@ void mpmont_destroy(mpmont *mm) MP_DROP(mm->m); MP_DROP(mm->r); MP_DROP(mm->r2); + MP_DROP(mm->mi); } /* --- @mpmont_reduce@ --- * * * 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) { - mpw *dv, *dvl; - const mpw *mv, *mvl; - size_t n; + mp_div(0, &d, a, mm->m); + return (d); +} + +#else - /* --- Initial conditioning of the arguments --- */ +mp *mpmont_reduce(mpmont *mm, mp *d, mp *a) +{ + size_t n = mm->n; + + /* --- Check for serious Karatsuba reduction --- */ + + if (n > MPK_THRESH * 3) { + mp al; + mpw *vl; + mp *u; + + if (MP_LEN(a) >= n) + vl = a->v + n; + else + vl = a->vl; + mp_build(&al, a->v, vl); + u = mp_mul(MP_NEW, &al, mm->mi); + if (MP_LEN(u) > n) + u->vl = u->v + n; + u = mp_mul(u, u, mm->m); + d = mp_add(d, a, u); + mp_drop(u); + } - n = MP_LEN(mm->m); + /* --- Otherwise do it the hard way --- */ - if (d == a) - MP_MODIFY(d, 2 * n); else { - MP_MODIFY(d, 2 * n); - memcpy(d->v, a->v, MPWS(MP_LEN(a))); - memset(d->v + MP_LEN(a), 0, MPWS(MP_LEN(d) - MP_LEN(a))); - } - - dv = d->v; dvl = d->vl; - mv = mm->m->v; mvl = mm->m->vl; + mpw *dv, *dvl; + mpw *mv, *mvl; + mpw mi; + size_t k = n; + + /* --- Initial conditioning of the arguments --- */ + + a = MP_COPY(a); + if (d) + MP_DROP(d); + d = a; + MP_DEST(d, 2 * n + 1, a->f); - /* --- Let's go to work --- */ + dv = d->v; dvl = d->vl; + mv = mm->m->v; mvl = mm->m->vl; - while (n--) { - mpw u = MPW(*dv * mm->mi); - MPX_UMLAN(dv, dvl, mv, mvl, u); - dv++; + /* --- Let's go to work --- */ + + mi = mm->mi->v[0]; + while (k--) { + mpw u = MPW(*dv * mi); + MPX_UMLAN(dv, dvl, mv, mvl, u); + dv++; + } } - /* --- Done --- */ + /* --- Wrap everything up --- */ - memmove(d->v, dv, MPWS(dvl - dv)); - d->vl -= dv - d->v; + memmove(d->v, d->v + n, MPWS(MP_LEN(d) - n)); + d->vl -= n; + if (MPX_UCMP(d->v, d->vl, >=, mm->m->v, mm->m->vl)) + mpx_usub(d->v, d->vl, d->v, d->vl, mm->m->v, mm->m->vl); + if (d->f & MP_NEG) { + mpx_usub(d->v, d->vl, mm->m->v, mm->m->vl, d->v, d->vl); + d->f &= ~MP_NEG; + } MP_SHRINK(d); - d->f = a->f & MP_BURN; - if (MP_CMP(d, >=, mm->m)) - d = mp_sub(d, d, mm->m); 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; - - /* --- 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++; - } +mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b) +{ + if (mm->n > MPK_THRESH * 3) { + 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; + mpw mi; - /* --- Simpler Montgomery reduction phase --- */ + /* --- Initial conditioning of the arguments --- */ - while (i < n) { - mpw u = MPW(*dv * mm->mi); - MPX_UMLAN(dv, dvl, mv, mvl, u); - dv++; - i++; - } + 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_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; + bv = b->v; bvl = b->vl; + mv = mm->m->v; mvl = mm->m->vl; + y = *bv; + + /* --- Montgomery multiplication phase --- */ + + i = 0; + mi = mm->mi->v[0]; + while (i < n && av < avl) { + mpw x = *av++; + mpw u = MPW((*dv + x * y) * mi); + MPX_UMLAN(dv, dvl, bv, bvl, x); + MPX_UMLAN(dv, dvl, mv, mvl, u); + dv++; + i++; + } - /* --- Done --- */ + /* --- Simpler Montgomery reduction phase --- */ + + while (i < n) { + mpw u = MPW(*dv * mi); + MPX_UMLAN(dv, dvl, mv, mvl, u); + dv++; + i++; + } + + /* --- Done --- */ + + memmove(d->v, dv, MPWS(dvl - dv)); + d->vl -= dv - d->v; + if (MPX_UCMP(d->v, d->vl, >=, mm->m->v, mm->m->vl)) + mpx_usub(d->v, d->vl, d->v, d->vl, mm->m->v, mm->m->vl); + if ((a->f ^ b->f) & MP_NEG) + mpx_usub(d->v, d->vl, mm->m->v, mm->m->vl, d->v, d->vl); + MP_SHRINK(d); + d->f = (a->f | b->f) & MP_BURN; + MP_DROP(a); + MP_DROP(b); + } - 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 + +/*----- Exponentiation ----------------------------------------------------*/ + /* --- @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$%. + * Returns: Result, %$(a R^{-1})^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 *spare = MP_NEW; - - mp_scan(&sc, e); - - if (MP_STEP(&sc)) { - size_t sq = 0; - for (;;) { - mp *dd; - if (MP_BIT(&sc)) { - while (sq) { - dd = mpmont_mul(mm, spare, ar, ar); - spare = ar; ar = dd; - sq--; - } - dd = mpmont_mul(mm, spare, d, ar); - spare = d; d = dd; - } - sq++; - if (!MP_STEP(&sc)) - break; - } - } - MP_DROP(ar); - if (spare != MP_NEW) - MP_DROP(spare); - return (d); + mp *x = MP_COPY(mm->r); + mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW; + + MP_SHRINK(e); + if (MP_LEN(e) == 0) + ; + else if (MP_LEN(e) < EXP_THRESH) + EXP_SIMPLE(x, a, e); + else + EXP_WINDOW(x, a, e); + mp_drop(d); + mp_drop(spare); + 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_mul(mm, d, a, mm->r2); + d = mpmont_expr(mm, d, d, e); d = mpmont_reduce(mm, d, d); return (d); } @@ -306,15 +430,15 @@ static int tcreate(dstr *v) mpmont_create(&mm, m); - if (mm.mi != mi->v[0]) { + if (mm.mi->v[0] != mi->v[0]) { fprintf(stderr, "\n*** bad mi: found %lu, expected %lu", - (unsigned long)mm.mi, (unsigned long)mi->v[0]); + (unsigned long)mm.mi->v[0], (unsigned long)mi->v[0]); fputs("\nm = ", stderr); mp_writefile(m, stderr, 10); fputc('\n', stderr); ok = 0; } - if (MP_CMP(mm.r, !=, r)) { + if (!MP_EQ(mm.r, r)) { fputs("\n*** bad r", stderr); fputs("\nm = ", stderr); mp_writefile(m, stderr, 10); fputs("\nexpected ", stderr); mp_writefile(r, stderr, 10); @@ -323,7 +447,7 @@ static int tcreate(dstr *v) ok = 0; } - if (MP_CMP(mm.r2, !=, r2)) { + if (!MP_EQ(mm.r2, r2)) { fputs("\n*** bad r2", stderr); fputs("\nm = ", stderr); mp_writefile(m, stderr, 10); fputs("\nexpected ", stderr); mp_writefile(r2, stderr, 10); @@ -337,6 +461,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 +471,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_EQ(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_EQ(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,9 +535,9 @@ 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)) { + if (!MP_EQ(mr, r)) { fputs("\n*** montgomery modexp failed", stderr); fputs("\n m = ", stderr); mp_writefile(m, stderr, 10); fputs("\n a = ", stderr); mp_writefile(a, stderr, 10); @@ -426,14 +554,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 } }, };