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

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