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

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