From: mdw Date: Fri, 19 Nov 1999 13:19:29 +0000 (+0000) Subject: Simultaneous exponentiation support. X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/commitdiff_plain/4cc02d06a57cc386ba445d34d9b83329bb29f89c Simultaneous exponentiation support. --- diff --git a/mpmont-mexp.c b/mpmont-mexp.c new file mode 100644 index 0000000..4180fbc --- /dev/null +++ b/mpmont-mexp.c @@ -0,0 +1,299 @@ +/* -*-c-*- + * + * $Id: mpmont-mexp.c,v 1.1 1999/11/19 13:19:29 mdw Exp $ + * + * Multiplle simultaneous exponentiations + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mpmont-mexp.c,v $ + * Revision 1.1 1999/11/19 13:19:29 mdw + * Simultaneous exponentiation support. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "mp.h" +#include "mpmont.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @mpmont_mexpr@ --- * + * + * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context + * @mpmont_factor *f@ = pointer to array of factors + * @size_t n@ = number of factors supplied + * + * Returns: If the bases are %$g_0, g_1, \ldots, g_{n-1}$% and the + * exponents are %$e_0, e_1, \ldots, e_{n-1}$% then the result + * is: + * + * %$g_0^{e_0} g_1^{e_1} \ldots g_{n-1}^{e_{n-1}} R \bmod m$% + */ + +typedef struct scan { + size_t len; + mpw w; +} scan; + +mp *mpmont_mexpr(mpmont *mm, mpmont_factor *f, size_t n) +{ + size_t vn = 1 << n; + mp **v = xmalloc(vn * sizeof(mp *)); + scan *s; + size_t o; + unsigned b; + mp *a = MP_COPY(mm->r); + mp *spare = MP_NEW; + + /* --- Perform the precomputation --- */ + + { + size_t i, j; + size_t mask; + + /* --- Fill in the rest of the array --- * + * + * Zero never gets used. + */ + + j = 0; + mask = 0; + for (i = 1; i < vn; i++) { + + /* --- Check for a new bit entering --- * + * + * If a bit gets set that wasn't set before, then all the lower bits + * are zeroes and I've got to introduce a new base into the array. + */ + + if ((i & mask) == 0) { + v[i] = mpmont_mul(mm, MP_NEW, f[j++].base, mm->r2); + mask = i; + } + + /* --- Otherwise I can get away with a single multiplication --- * + * + * In particular, if %$i$% has more than one bit set, then I only need + * to calculate %$v_i = v_{\mathit{mask}} v_{i - \mathit{mask}}$%. + * Since both are less than %$i$%, they must have already been + * computed. + */ + + else + v[i] = mpmont_mul(mm, MP_NEW, v[mask], v[i & ~mask]); + } + } + + /* --- Set up the bitscanners --- * + * + * I must scan the exponents from left to right, which is a shame. It + * means that I can't use the standard @mpscan@ stuff, in particular. + */ + + { + size_t i; + + s = xmalloc(n * sizeof(scan)); + o = 0; + for (i = 0; i < n; i++) { + s[i].len = MP_LEN(f[i].exp); + if (s[i].len > o) + o = s[i].len; + } + b = 0; + } + + /* --- Now do the actual calculation --- */ + + b = 0; + for (;;) { + size_t i; + size_t j; + mp *dd; + + /* --- If no more bits, get some more --- */ + + if (!b) { + if (!o) + break; + o--; + b = MPW_BITS; + } + + /* --- Work out the next index --- */ + + j = 0; + b--; + for (i = 0; i < n; i++) { + if (o < s[i].len) + j |= (((f[i].exp->v[o] >> b) & 1) << i); + } + + /* --- Accumulate the result --- */ + + if (spare) { + dd = mpmont_mul(mm, spare, a, a); + spare = a; + a = dd; + } + + if (j) { + dd = mpmont_mul(mm, spare, a, v[j]); + spare = a; + a = dd; + } + } + + /* --- Tidy up afterwards --- */ + + { + size_t i; + for (i = 1; i < vn; i++) + MP_DROP(v[i]); + if (spare) + MP_DROP(spare); + free(v); + free(s); + } + + return (a); +} + +/* --- @mpmont_mexp@ --- * + * + * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context + * @mpmont_factor *f@ = pointer to array of factors + * @size_t n@ = number of factors supplied + * + * Returns: Product of bases raised to exponents, all mod @m@. + * + * Use: Convenient interface over @mpmont_mexpr@. + */ + +mp *mpmont_mexp(mpmont *mm, mpmont_factor *f, size_t n) +{ + mp *d = mpmont_mexpr(mm, f, n); + d = mpmont_reduce(mm, d, d); + return (d); +} + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +#include + +static int verify(size_t n, dstr *v) +{ + mp *m = *(mp **)v[0].buf; + mpmont_factor *f = xmalloc(n * sizeof(*f)); + mp *r, *rr; + size_t i, j; + mpmont mm; + int ok = 1; + + j = 1; + for (i = 0; i < n; i++) { + f[i].base = *(mp **)v[j++].buf; + f[i].exp = *(mp **)v[j++].buf; + } + + rr = *(mp **)v[j].buf; + mpmont_create(&mm, m); + r = mpmont_mexp(&mm, f, n); + if (MP_CMP(r, !=, rr)) { + fputs("\n*** mexp failed\n", stderr); + fputs("m = ", stderr); mp_writefile(m, stderr, 10); + for (i = 0; i < n; i++) { + fprintf(stderr, "\ng_%u = ", i); + mp_writefile(f[i].base, stderr, 10); + fprintf(stderr, "\ne_%u = ", i); + mp_writefile(f[i].exp, stderr, 10); + } + fputs("\nr = ", stderr); mp_writefile(r, stderr, 10); + fputs("\nR = ", stderr); mp_writefile(rr, stderr, 10); + fputc('\n', stderr); + ok = 0; + } + + for (i = 0; i < n; i++) { + MP_DROP(f[i].base); + MP_DROP(f[i].exp); + } + MP_DROP(m); + MP_DROP(r); + MP_DROP(rr); + mpmont_destroy(&mm); + return (ok); +} + +static int t1(dstr *v) { return verify(1, v); } +static int t2(dstr *v) { return verify(2, v); } +static int t3(dstr *v) { return verify(3, v); } +static int t4(dstr *v) { return verify(4, v); } +static int t5(dstr *v) { return verify(5, v); } + +static test_chunk tests[] = { + { "mexp-1", t1, { &type_mp, + &type_mp, &type_mp, + &type_mp, 0 } }, + { "mexp-2", t2, { &type_mp, + &type_mp, &type_mp, + &type_mp, &type_mp, + &type_mp, 0 } }, + { "mexp-3", t3, { &type_mp, + &type_mp, &type_mp, + &type_mp, &type_mp, + &type_mp, &type_mp, + &type_mp, 0 } }, + { "mexp-4", t4, { &type_mp, + &type_mp, &type_mp, + &type_mp, &type_mp, + &type_mp, &type_mp, + &type_mp, &type_mp, + &type_mp, 0 } }, + { "mexp-5", t5, { &type_mp, + &type_mp, &type_mp, + &type_mp, &type_mp, + &type_mp, &type_mp, + &type_mp, &type_mp, + &type_mp, &type_mp, + &type_mp, 0 } }, + { 0, 0, { 0 } } +}; + +int main(int argc, char *argv[]) +{ + sub_init(); + test_run(argc, argv, tests, SRCDIR "/tests/mpmont"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/