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