chiark / gitweb /
Merge branch '2.4.x' into 2.5.x
[catacomb] / 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