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