git.distorted.org.uk Git - u/mdw/catacomb/blob - mpmont-mexp.c

   1 /* -*-c-*-
   2  *
   3  * $Id: mpmont-mexp.c,v 1.2 1999/11/21 11:35:10 mdw Exp $
   4  *
   5  * Multiplle simultaneous exponentiations
   6  *
   7  * (c) 1999 Straylight/Edgeware
   8  */
   9
  10 /*----- Licensing notice --------------------------------------------------*
  11  *
  12  * This file is part of Catacomb.
  13  *
  14  * Catacomb is free software; you can redistribute it and/or modify
  15  * it under the terms of the GNU Library General Public License as
  16  * published by the Free Software Foundation; either version 2 of the
  17  * License, or (at your option) any later version.
  18  *
  19  * Catacomb is distributed in the hope that it will be useful,
  20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22  * GNU Library General Public License for more details.
  23  *
  24  * You should have received a copy of the GNU Library General Public
  25  * License along with Catacomb; if not, write to the Free
  26  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  27  * MA 02111-1307, USA.
  28  */
  29
  30 /*----- Revision history --------------------------------------------------*
  31  *
  32  * $Log: mpmont-mexp.c,v $
  33  * Revision 1.2  1999/11/21 11:35:10  mdw
  34  * Performance improvement: use @mp_sqr@ and @mpmont_reduce@ instead of
  35  * @mpmont_mul@ for squaring in exponentiation.
  36  *
  37  * Revision 1.1  1999/11/19 13:19:29  mdw
  38  * Simultaneous exponentiation support.
  39  *
  40  */
  41
  42 /*----- Header files ------------------------------------------------------*/
  43
  44 #include "mp.h"
  45 #include "mpmont.h"
  46
  47 /*----- Main code ---------------------------------------------------------*/
  48
  49 /* --- @mpmont_mexpr@ --- *
  50  *
  51  * Arguments:   @mpmont *mm@ = pointer to Montgomery reduction context
  52  *              @mpmont_factor *f@ = pointer to array of factors
  53  *              @size_t n@ = number of factors supplied
  54  *
  55  * Returns:     If the bases are %$g_0, g_1, \ldots, g_{n-1}$% and the
  56  *              exponents are %$e_0, e_1, \ldots, e_{n-1}$% then the result
  57  *              is:
  58  *
  59  *              %$g_0^{e_0} g_1^{e_1} \ldots g_{n-1}^{e_{n-1}} R \bmod m$%
  60  */
  61
  62 typedef struct scan {
  63   size_t len;
  64   mpw w;
  65 } scan;
  66
  67 mp *mpmont_mexpr(mpmont *mm, mpmont_factor *f, size_t n)
  68 {
  69   size_t vn = 1 << n;
  70   mp **v = xmalloc(vn * sizeof(mp *));
  71   scan *s;
  72   size_t o;
  73   unsigned b;
  74   mp *a = MP_COPY(mm->r);
  75   mp *spare = MP_NEW;
  76
  77   /* --- Perform the precomputation --- */
  78
  79   {
  80     size_t i, j;
  81     size_t mask;
  82
  83     /* --- Fill in the rest of the array --- *
  84      *
  85      * Zero never gets used.
  86      */
  87
  88     j = 0;
  89     mask = 0;
  90     for (i = 1; i < vn; i++) {
  91
  92       /* --- Check for a new bit entering --- *
  93        *
  94        * If a bit gets set that wasn't set before, then all the lower bits
  95        * are zeroes and I've got to introduce a new base into the array.
  96        */
  97
  98       if ((i & mask) == 0) {
  99         v[i] = mpmont_mul(mm, MP_NEW, f[j++].base, mm->r2);
 100         mask = i;
 101       }
 102
 103       /* --- Otherwise I can get away with a single multiplication --- *
 104        *
 105        * In particular, if %$i$% has more than one bit set, then I only need
 106        * to calculate %$v_i = v_{\mathit{mask}} v_{i - \mathit{mask}}$%.
 107        * Since both are less than %$i$%, they must have already been
 108        * computed.
 109        */
 110
 111       else
 112         v[i] = mpmont_mul(mm, MP_NEW, v[mask], v[i & ~mask]);
 113     }
 114   }
 115
 116   /* --- Set up the bitscanners --- *
 117    *
 118    * I must scan the exponents from left to right, which is a shame.  It
 119    * means that I can't use the standard @mpscan@ stuff, in particular.
 120    */
 121
 122   {
 123     size_t i;
 124
 125     s = xmalloc(n * sizeof(scan));
 126     o = 0;
 127     for (i = 0; i < n; i++) {
 128       s[i].len = MP_LEN(f[i].exp);
 129       if (s[i].len > o)
 130         o = s[i].len;
 131     }
 132     b = 0;
 133   }
 134
 135   /* --- Now do the actual calculation --- */
 136
 137   b = 0;
 138   for (;;) {
 139     size_t i;
 140     size_t j;
 141     mp *dd;
 142
 143     /* --- If no more bits, get some more --- */
 144
 145     if (!b) {
 146       if (!o)
 147         break;
 148       o--;
 149       b = MPW_BITS;
 150     }
 151
 152     /* --- Work out the next index --- */
 153
 154     j = 0;
 155     b--;
 156     for (i = 0; i < n; i++) {
 157       if (o < s[i].len)
 158         j |= (((f[i].exp->v[o] >> b) & 1) << i);
 159     }
 160
 161     /* --- Accumulate the result --- */
 162
 163     if (spare) {
 164       dd = mp_sqr(spare, a);
 165       dd = mpmont_reduce(mm, dd, dd);
 166       spare = a;
 167       a = dd;
 168     }
 169
 170     if (j) {
 171       dd = mpmont_mul(mm, spare, a, v[j]);
 172       spare = a;
 173       a = dd;
 174     }
 175   }
 176
 177   /* --- Tidy up afterwards --- */
 178
 179   {
 180     size_t i;
 181     for (i = 1; i < vn; i++)
 182       MP_DROP(v[i]);
 183     if (spare)
 184       MP_DROP(spare);
 185     free(v);
 186     free(s);
 187   }
 188
 189   return (a);
 190 }
 191
 192 /* --- @mpmont_mexp@ --- *
 193  *
 194  * Arguments:   @mpmont *mm@ = pointer to Montgomery reduction context
 195  *              @mpmont_factor *f@ = pointer to array of factors
 196  *              @size_t n@ = number of factors supplied
 197  *
 198  * Returns:     Product of bases raised to exponents, all mod @m@.
 199  *
 200  * Use:         Convenient interface over @mpmont_mexpr@.
 201  */
 202
 203 mp *mpmont_mexp(mpmont *mm, mpmont_factor *f, size_t n)
 204 {
 205   mp *d = mpmont_mexpr(mm, f, n);
 206   d = mpmont_reduce(mm, d, d);
 207   return (d);
 208 }
 209
 210 /*----- Test rig ----------------------------------------------------------*/
 211
 212 #ifdef TEST_RIG
 213
 214 #include <mLib/testrig.h>
 215
 216 static int verify(size_t n, dstr *v)
 217 {
 218   mp *m = *(mp **)v[0].buf;
 219   mpmont_factor *f = xmalloc(n * sizeof(*f));
 220   mp *r, *rr;
 221   size_t i, j;
 222   mpmont mm;
 223   int ok = 1;
 224
 225   j = 1;
 226   for (i = 0; i < n; i++) {
 227     f[i].base = *(mp **)v[j++].buf;
 228     f[i].exp = *(mp **)v[j++].buf;
 229   }
 230
 231   rr = *(mp **)v[j].buf;
 232   mpmont_create(&mm, m);
 233   r = mpmont_mexp(&mm, f, n);
 234   if (MP_CMP(r, !=, rr)) {
 235     fputs("\n*** mexp failed\n", stderr);
 236     fputs("m = ", stderr); mp_writefile(m, stderr, 10);
 237     for (i = 0; i < n; i++) {
 238       fprintf(stderr, "\ng_%u = ", i);
 239       mp_writefile(f[i].base, stderr, 10);
 240       fprintf(stderr, "\ne_%u = ", i);
 241       mp_writefile(f[i].exp, stderr, 10);
 242     }
 243     fputs("\nr = ", stderr); mp_writefile(r, stderr, 10);
 244     fputs("\nR = ", stderr); mp_writefile(rr, stderr, 10);
 245     fputc('\n', stderr);
 246     ok = 0;
 247   }
 248
 249   for (i = 0; i < n; i++) {
 250     MP_DROP(f[i].base);
 251     MP_DROP(f[i].exp);
 252   }
 253   MP_DROP(m);
 254   MP_DROP(r);
 255   MP_DROP(rr);
 256   mpmont_destroy(&mm);
 257   return (ok);
 258 }
 259
 260 static int t1(dstr *v) { return verify(1, v); }
 261 static int t2(dstr *v) { return verify(2, v); }
 262 static int t3(dstr *v) { return verify(3, v); }
 263 static int t4(dstr *v) { return verify(4, v); }
 264 static int t5(dstr *v) { return verify(5, v); }
 265
 266 static test_chunk tests[] = {
 267   { "mexp-1", t1, { &type_mp,
 268                     &type_mp, &type_mp,
 269                     &type_mp, 0 } },
 270   { "mexp-2", t2, { &type_mp,
 271                     &type_mp, &type_mp,
 272                     &type_mp, &type_mp,
 273                     &type_mp, 0 } },
 274   { "mexp-3", t3, { &type_mp,
 275                     &type_mp, &type_mp,
 276                     &type_mp, &type_mp,
 277                     &type_mp, &type_mp,
 278                     &type_mp, 0 } },
 279   { "mexp-4", t4, { &type_mp,
 280                     &type_mp, &type_mp,
 281                     &type_mp, &type_mp,
 282                     &type_mp, &type_mp,
 283                     &type_mp, &type_mp,
 284                     &type_mp, 0 } },
 285   { "mexp-5", t5, { &type_mp,
 286                     &type_mp, &type_mp,
 287                     &type_mp, &type_mp,
 288                     &type_mp, &type_mp,
 289                     &type_mp, &type_mp,
 290                     &type_mp, &type_mp,
 291                     &type_mp, 0 } },
 292   { 0, 0, { 0 } }
 293 };
 294
 295 int main(int argc, char *argv[])
 296 {
 297   sub_init();
 298   test_run(argc, argv, tests, SRCDIR "/tests/mpmont");
 299   return (0);
 300 }
 301
 302 #endif
 303
 304 /*----- That's all, folks -------------------------------------------------*/