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

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