chiark - git - mdw - catacomb/blob - mpmont-mexp.c

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