chiark - git - mdw - catacomb/blob - exp.h

   1 /* -*-c-*-
   2  *
   3  * $Id$
   4  *
   5  * Generalized exponentiation
   6  *
   7  * (c) 2001 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 #ifdef CATACOMB_EXP_H
  31 #  error "Multiple inclusion of <catacomb/exp.h>"
  32 #endif
  33
  34 #define CATACOMB_EXP_H
  35
  36 #ifdef __cplusplus
  37   extern "C" {
  38 #endif
  39
  40 /*----- Header files ------------------------------------------------------*/
  41
  42 #include <stddef.h>
  43
  44 #include <mLib/alloc.h>
  45
  46 #ifndef CATACOMB_MP_H
  47 #  include "mp.h"
  48 #endif
  49
  50 /*----- Data structures ---------------------------------------------------*/
  51
  52 typedef struct exp_simulscan {
  53   mpw w;
  54   size_t len;
  55   const mpw *v;
  56 } exp_simulscan;
  57
  58 typedef struct exp_simul {
  59   unsigned b;
  60   size_t o, n;
  61   exp_simulscan *s;
  62 } exp_simul;
  63
  64 /*----- Macros provided ---------------------------------------------------*/
  65
  66 /* --- Parameters --- */
  67
  68 #ifndef EXP_WINSZ                       /* Sliding window size */
  69 #  define EXP_WINSZ 4                   /* Predefine if you need to */
  70 #endif
  71
  72 /* --- These are determined from the window size --- *
  73  *
  74  * Given a %$k$%-bit exponent, I expect to do %$k/2$% multiplies if I use the
  75  * simple way.  If I use an n-bit sliding window, then I do %$2^n$%
  76  * multiplies up front, but I only do %$(2^n - 1)/2^n k/n$% multiplies for
  77  * the exponentiation.  This is a win when
  78  *
  79  *   %$k \ge \frac{n 2^{n+1}}{n - 2}$%
  80  */
  81
  82 #define EXP_TABSZ (1 << EXP_WINSZ)
  83 #define EXP_THRESH                                                      \
  84     ((EXP_WINSZ * (2 << EXP_WINSZ))/((EXP_WINSZ - 2) * MPW_BITS))
  85
  86 /* --- Required operations --- *
  87  *
  88  * The macros here are independent of the underlying group elements.  You
  89  * must provide the necessary group operations and other definitions.  The
  90  * group operation is assumed to be written multiplicatively.
  91  *
  92  * @EXP_TYPE@                   The type of a group element, e.g., @mp *@.
  93  *
  94  * @EXP_COPY(d, x)@             Makes @d@ be a copy of @x@.
  95  *
  96  * @EXP_DROP(x)@                Discards the element @x@, reclaiming any
  97  *                              memory it used.
  98  *
  99  * @EXP_MUL(a, x)@              Multiplies @a@ by @x@ (writing the result
 100  *                              back to @a@).
 101  *
 102  * @EXP_FIX(x)@                 Makes @x@ be a canonical representation of
 103  *                              its value.  All multiplications have the
 104  *                              right argument canonical.
 105  *
 106  * @EXP_SQR(a)@                 Multiplies @a@ by itself.
 107  *
 108  * @EXP_SETMUL(d, x, y)@        Sets @d@ to be the product of @x@ and @y@.
 109  *                              The value @d@ has not been initialized.
 110  *
 111  * @EXP_SETSQR(d, x)@           Sets @d@ to be the square of @x@.
 112  *
 113  * Only @EXP_TYPE@, @EXP_MUL@ and @EXP_SQR@ are required for simple
 114  * exponentation.  Sliding window and simultaneous exponentation require all
 115  * of the operations.
 116  */
 117
 118 #ifndef EXP_TYPE
 119 #  error "EXP_TYPE not defined for <catacomb/exp.h>"
 120 #endif
 121
 122 /* --- @EXP_SIMPLE@ --- *
 123  *
 124  * Arguments:   @a@ = the result object, initially a multiplicative identity
 125  *              @g@ = the object to exponentiate
 126  *              @x@ = the exponent, as a multiprecision integer
 127  *
 128  * Use:         Performs a simple left-to-right exponentiation.  At the end
 129  *              of the code, the answer is left in @a@; @g@ and @x@ are
 130  *              unchanged.
 131  */
 132
 133 #define EXP_SIMPLE(a, g, x) do {                                        \
 134   mpscan sc;                                                            \
 135   unsigned sq = 0;                                                      \
 136                                                                         \
 137   /* --- Begin scanning --- */                                          \
 138                                                                         \
 139   mp_rscan(&sc, x);                                                     \
 140   if (!MP_RSTEP(&sc))                                                   \
 141     goto exp_simple_exit;                                               \
 142   while (!MP_RBIT(&sc))                                                 \
 143     MP_RSTEP(&sc);                                                      \
 144                                                                         \
 145   /* --- Do the main body of the work --- */                            \
 146                                                                         \
 147   EXP_FIX(g);                                                           \
 148   for (;;) {                                                            \
 149     EXP_MUL(a, g);                                                      \
 150     sq = 0;                                                             \
 151     for (;;) {                                                          \
 152       if (!MP_RSTEP(&sc))                                               \
 153         goto exp_simple_done;                                           \
 154       sq++;                                                             \
 155       if (MP_RBIT(&sc))                                                 \
 156         break;                                                          \
 157     }                                                                   \
 158     while (sq--) EXP_SQR(a);                                            \
 159   }                                                                     \
 160                                                                         \
 161   /* --- Do a final round of squaring --- */                            \
 162                                                                         \
 163 exp_simple_done:                                                        \
 164   while (sq--) EXP_SQR(a);                                              \
 165 exp_simple_exit:;                                                       \
 166 } while (0)
 167
 168 /* --- @EXP_WINDOW@ --- *
 169  *
 170  * Arguments:   @a@ = the result object, initially a multiplicative identity
 171  *              @g@ = the object to exponentiate
 172  *              @x@ = the exponent, as a multiprecision integer
 173  *
 174  * Use:         Performs a sliding-window exponentiation.  At the end of the
 175  *              code, the answer is left in @a@; @g@ and @x@ are unchanged.
 176  */
 177
 178 #define EXP_WINDOW(a, g, x) do {                                        \
 179   EXP_TYPE *v;                                                          \
 180   EXP_TYPE g2;                                                          \
 181   unsigned i, sq = 0;                                                   \
 182   mpscan sc;                                                            \
 183                                                                         \
 184   /* --- Get going --- */                                               \
 185                                                                         \
 186   mp_rscan(&sc, x);                                                     \
 187   if (!MP_RSTEP(&sc))                                                   \
 188     goto exp_window_exit;                                               \
 189                                                                         \
 190   /* --- Do the precomputation --- */                                   \
 191                                                                         \
 192   EXP_FIX(g);                                                           \
 193   EXP_SETSQR(g2, g);                                                    \
 194   EXP_FIX(g2);                                                          \
 195   v = xmalloc(EXP_TABSZ * sizeof(EXP_TYPE));                            \
 196   EXP_COPY(v[0], g);                                                    \
 197   for (i = 1; i < EXP_TABSZ; i++) {                                     \
 198     EXP_SETMUL(v[i], v[i - 1], g2);                                     \
 199     EXP_FIX(v[i]);                                                      \
 200   }                                                                     \
 201   EXP_DROP(g2);                                                         \
 202                                                                         \
 203   /* --- Skip top-end zero bits --- *                                   \
 204    *                                                                    \
 205    * If the initial step worked, there must be a set bit somewhere, so  \
 206    * keep stepping until I find it.                                     \
 207    */                                                                   \
 208                                                                         \
 209   while (!MP_RBIT(&sc))                                                 \
 210     MP_RSTEP(&sc);                                                      \
 211                                                                         \
 212   /* --- Now for the main work --- */                                   \
 213                                                                         \
 214   for (;;) {                                                            \
 215     unsigned l = 1;                                                     \
 216     unsigned z = 0;                                                     \
 217                                                                         \
 218     /* --- The next bit is set, so read a window index --- *            \
 219      *                                                                  \
 220      * Reset @i@ to zero and increment @sq@.  Then, until either I read \
 221      * @WINSZ@ bits or I run out of bits, scan in a bit: if it's clear, \
 222      * bump the @z@ counter; if it's set, push a set bit into @i@,      \
 223      * shift it over by @z@ bits, bump @sq@ by @z + 1@ and clear @z@.   \
 224      * By the end of this palaver, @i@ is an index to the precomputed   \
 225      * value in @v@.                                                    \
 226      */                                                                 \
 227                                                                         \
 228     i = 0;                                                              \
 229     sq++;                                                               \
 230     while (l < EXP_WINSZ && MP_RSTEP(&sc)) {                            \
 231       l++;                                                              \
 232       if (!MP_RBIT(&sc))                                                \
 233         z++;                                                            \
 234       else {                                                            \
 235         i = ((i << 1) | 1) << z;                                        \
 236         sq += z + 1;                                                    \
 237         z = 0;                                                          \
 238       }                                                                 \
 239     }                                                                   \
 240                                                                         \
 241     /* --- Do the squaring --- *                                        \
 242      *                                                                  \
 243      * Remember that @sq@ carries over from the zero-skipping stuff     \
 244      * below.                                                           \
 245      */                                                                 \
 246                                                                         \
 247     while (sq--) EXP_SQR(a);                                            \
 248                                                                         \
 249     /* --- Do the multiply --- */                                       \
 250                                                                         \
 251     EXP_MUL(a, v[i]);                                                   \
 252                                                                         \
 253     /* --- Now grind along through the rest of the bits --- */          \
 254                                                                         \
 255     sq = z;                                                             \
 256     for (;;) {                                                          \
 257       if (!MP_RSTEP(&sc))                                               \
 258         goto exp_window_done;                                           \
 259       if (MP_RBIT(&sc))                                                 \
 260         break;                                                          \
 261       sq++;                                                             \
 262     }                                                                   \
 263   }                                                                     \
 264                                                                         \
 265   /* --- Do a final round of squaring --- */                            \
 266                                                                         \
 267 exp_window_done:                                                        \
 268   while (sq--) EXP_SQR(a);                                              \
 269   for (i = 0; i < EXP_TABSZ; i++)                                       \
 270     EXP_DROP(v[i]);                                                     \
 271   xfree(v);                                                             \
 272 exp_window_exit:;                                                       \
 273 } while (0)
 274
 275 /* --- @EXP_SIMUL@ --- *
 276  *
 277  * Arguments:   @a@ = the result object, initially a multiplicative identity
 278  *              @f@ = pointer to a vector of base/exp pairs
 279  *              @n@ = the number of base/exp pairs
 280  *
 281  * Use:         Performs a simultaneous sliding-window exponentiation.  The
 282  *              @f@ table is an array of structures containing members @base@
 283  *              of type @EXP_TYPE@, and @exp@ of type @mp *@.
 284  */
 285
 286 #define EXP_SIMUL(a, f, n) do {                                         \
 287   size_t i, j, jj, k;                                                   \
 288   size_t vn = 1 << (EXP_WINSZ * n), m = (1 << n) - 1;                   \
 289   EXP_TYPE *v = xmalloc(vn * sizeof(EXP_TYPE));                         \
 290   exp_simul e;                                                          \
 291   unsigned sq = 0;                                                      \
 292                                                                         \
 293   /* --- Fill in the precomputed table --- */                           \
 294                                                                         \
 295   j = 1;                                                                \
 296   for (i = 0; i < n; i++) {                                             \
 297     EXP_COPY(v[j], f[n - 1 - i].base);                                  \
 298     EXP_FIX(v[j]);                                                      \
 299     j <<= 1;                                                            \
 300   }                                                                     \
 301   k = n * EXP_WINSZ;                                                    \
 302   jj = 1;                                                               \
 303   for (; i < k; i++) {                                                  \
 304     EXP_SETSQR(v[j], v[jj]);                                            \
 305     EXP_FIX(v[j]);                                                      \
 306     j <<= 1; jj <<= 1;                                                  \
 307   }                                                                     \
 308   for (i = 1; i < vn; i <<= 1) {                                        \
 309     for (j = 1; j < i; j++) {                                           \
 310       EXP_SETMUL(v[j + i], v[j], v[i]);                                 \
 311       EXP_FIX(v[j + i]);                                                \
 312     }                                                                   \
 313   }                                                                     \
 314                                                                         \
 315   /* --- Set up the bitscanners --- *                                   \
 316    *                                                                    \
 317    * Got to use custom scanners, to keep them all in sync.              \
 318    */                                                                   \
 319                                                                         \
 320   e.n = n;                                                              \
 321   e.b = 0;                                                              \
 322   e.s = xmalloc(n * sizeof(*e.s));                                      \
 323   e.o = 0;                                                              \
 324   for (i = 0; i < n; i++) {                                             \
 325     MP_SHRINK(f[i].exp);                                                \
 326     e.s[i].len = MP_LEN(f[i].exp);                                      \
 327     e.s[i].v = f[i].exp->v;                                             \
 328     if (e.s[i].len > e.o)                                               \
 329       e.o = e.s[i].len;                                                 \
 330   }                                                                     \
 331                                                                         \
 332   /* --- Skip as far as a nonzero column in the exponent matrix --- */  \
 333                                                                         \
 334   do {                                                                  \
 335     if (!e.o && !e.b)                                                   \
 336       goto exp_simul_done;                                              \
 337     i = exp_simulnext(&e, 0);                                           \
 338   } while (!(i & m));                                                   \
 339                                                                         \
 340   /* --- Now for the main work --- */                                   \
 341                                                                         \
 342   for (;;) {                                                            \
 343     unsigned l = 1;                                                     \
 344     unsigned z = 0;                                                     \
 345                                                                         \
 346     /* --- Just read a nonzero column, so read a window index --- *     \
 347      *                                                                  \
 348      * Clear high bits of @i@ and increment @sq@.  Then, until either I \
 349      * read @WINSZ@ columns or I run out, scan in a column and append   \
 350      * it to @i@.  If it's zero, bump the @z@ counter; if it's nonzero, \
 351      * bump @sq@ by @z + 1@ and clear @z@.  By the end of this palaver, \
 352      * @i@ is an index to the precomputed value in @v@, followed by     \
 353      * @n * z@ zero bits.                                               \
 354      */                                                                 \
 355                                                                         \
 356     sq++;                                                               \
 357     while (l < EXP_WINSZ && (e.o || e.b)) {                             \
 358       l++;                                                              \
 359       i = exp_simulnext(&e, i);                                         \
 360       if (!(i & m))                                                     \
 361         z++;                                                            \
 362       else {                                                            \
 363         sq += z + 1;                                                    \
 364         z = 0;                                                          \
 365       }                                                                 \
 366     }                                                                   \
 367                                                                         \
 368     /* --- Do the squaring --- *                                        \
 369      *                                                                  \
 370      * Remember that @sq@ carries over from the zero-skipping stuff     \
 371      * below.                                                           \
 372      */                                                                 \
 373                                                                         \
 374     while (sq--) EXP_SQR(a);                                            \
 375                                                                         \
 376     /* --- Do the multiply --- */                                       \
 377                                                                         \
 378     i >>= (z * n);                                                      \
 379     EXP_MUL(a, v[i]);                                                   \
 380                                                                         \
 381     /* --- Now grind along through the rest of the bits --- */          \
 382                                                                         \
 383     sq = z;                                                             \
 384     for (;;) {                                                          \
 385       if (!e.o && !e.b)                                                 \
 386         goto exp_simul_done;                                            \
 387       if ((i = exp_simulnext(&e, 0)) != 0)                              \
 388         break;                                                          \
 389       sq++;                                                             \
 390     }                                                                   \
 391   }                                                                     \
 392                                                                         \
 393   /* --- Do a final round of squaring --- */                            \
 394                                                                         \
 395 exp_simul_done:                                                         \
 396   while (sq--) EXP_SQR(a);                                              \
 397   for (i = 1; i < vn; i++)                                              \
 398     EXP_DROP(v[i]);                                                     \
 399   xfree(v);                                                             \
 400   xfree(e.s);                                                           \
 401 } while (0)
 402
 403 /*----- Functions provided ------------------------------------------------*/
 404
 405 /* --- @exp_simulnext@ --- *
 406  *
 407  * Arguments:   @exp_simul *e@ = pointer to state structure
 408  *              @size_t x@ = a current accumulator
 409  *
 410  * Returns:     The next column of bits.
 411  *
 412  * Use:         Scans the next column of bits for a simultaneous
 413  *              exponentiation.
 414  */
 415
 416 extern size_t exp_simulnext(exp_simul */*e*/, size_t /*x*/);
 417
 418 /*----- That's all, folks -------------------------------------------------*/
 419
 420 #ifdef __cplusplus
 421   }
 422 #endif