projects / u / mdw / catacomb / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
math/gfreduce.[ch]: Fix out-of-bounds memory access.
[u/mdw/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
Catacomb cryptographic library (downstream of http://git.distorted.org.uk/~mdw/catacomb/)