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