projects / u / mdw / catacomb / blob
? search:


7d31334f9da4796c775bbcf6d09779fb7d1f47b7
[u/mdw/catacomb] / mpreduce.c
1 /* -*-c-*-
2 *
3 * $Id: mpreduce.c,v 1.2 2004/04/08 01:36:15 mdw Exp $
4 *
5 * Efficient reduction modulo nice primes
6 *
7 * (c) 2004 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 /*----- Header files ------------------------------------------------------*/
31
32 #include <mLib/darray.h>
33 #include <mLib/macros.h>
34
35 #include "mp.h"
36 #include "mpreduce.h"
37 #include "mpreduce-exp.h"
38
39 /*----- Data structures ---------------------------------------------------*/
40
41 DA_DECL(instr_v, mpreduce_instr);
42
43 /*----- Main code ---------------------------------------------------------*/
44
45 /* --- @mpreduce_create@ --- *
46 *
47 * Arguments: @gfreduce *r@ = structure to fill in
48 * @mp *x@ = an integer
49 *
50 * Returns: ---
51 *
52 * Use: Initializes a context structure for reduction.
53 */
54
55 void mpreduce_create(mpreduce *r, mp *p)
56 {
57 mpscan sc;
58 enum { Z = 0, Z1 = 2, X = 4, X0 = 6 };
59 unsigned st = Z;
60 instr_v iv = DA_INIT;
61 unsigned long d, i;
62 unsigned op;
63 size_t w, b;
64
65 /* --- Fill in the easy stuff --- */
66
67 assert(MP_ISPOS(p));
68 d = mp_bits(p);
69 r->lim = d/MPW_BITS;
70 r->s = d%MPW_BITS;
71 if (r->s)
72 r->lim++;
73 r->p = mp_copy(p);
74
75 /* --- Stash a new instruction --- */
76
77 #define INSTR(op_, argx_, argy_) do { \
78 DA_ENSURE(&iv, 1); \
79 DA(&iv)[DA_LEN(&iv)].op = (op_); \
80 DA(&iv)[DA_LEN(&iv)].argx = (argx_); \
81 DA(&iv)[DA_LEN(&iv)].argy = (argy_); \
82 DA_EXTEND(&iv, 1); \
83 } while (0)
84
85 /* --- Main loop --- *
86 *
87 * A simple state machine decomposes @p@ conveniently into positive and
88 * negative powers of 2. The pure form of the state machine is left below
89 * for reference (and in case I need inspiration for a NAF exponentiator).
90 */
91
92 #ifdef DEBUG
93 for (i = 0, mp_scan(&sc, p); mp_step(&sc); i++) {
94 switch (st | mp_bit(&sc)) {
95 case Z | 1: st = Z1; break;
96 case Z1 | 0: st = Z; printf("+ %lu\n", i - 1); break;
97 case Z1 | 1: st = X; printf("- %lu\n", i - 1); break;
98 case X | 0: st = X0; break;
99 case X0 | 1: st = X; printf("- %lu\n", i - 1); break;
100 case X0 | 0: st = Z; printf("+ %lu\n", i - 1); break;
101 }
102 }
103 if (st >= X) printf("+ %lu\n", i);
104 #endif
105
106 for (i = 0, mp_scan(&sc, p); i < d - 1 && mp_step(&sc); i++) {
107 switch (st | mp_bit(&sc)) {
108 case Z | 1: st = Z1; break;
109 case Z1 | 0: st = Z; op = MPRI_SUB; goto instr;
110 case Z1 | 1: st = X; op = MPRI_ADD; goto instr;
111 case X | 0: st = X0; break;
112 case X0 | 1: st = X; op = MPRI_ADD; goto instr;
113 case X0 | 0: st = Z; op = MPRI_SUB; goto instr;
114 instr:
115 w = (d - i)/MPW_BITS + 1;
116 b = (MPW_BITS + i - d - 1)%MPW_BITS;
117 INSTR(op | !!b, w, b);
118 }
119 }
120 if ((DA(&iv)[DA_LEN(&iv) - 1].op & ~1u) == MPRI_SUB) {
121 fprintf(stderr,
122 "mpreduce can't handle moduli of the form 2^m + x\n");
123 abort();
124 }
125
126 #undef INSTR
127
128 /* --- Wrap up --- */
129
130 r->in = DA_LEN(&iv);
131 if (!r->s) {
132 r->iv = xmalloc(r->in * sizeof(mpreduce_instr));
133 memcpy(r->iv, DA(&iv), r->in * sizeof(mpreduce_instr));
134 } else {
135 r->iv = xmalloc(r->in * 2 * sizeof(mpreduce_instr));
136 for (i = 0; i < r->in; i++) {
137 r->iv[i] = DA(&iv)[i];
138 op = r->iv[i].op & ~1u;
139 w = r->iv[i].argx;
140 b = r->iv[i].argy;
141 b += r->s;
142 if (b >= MPW_BITS) {
143 b -= MPW_BITS;
144 w--;
145 }
146 if (b) op |= 1;
147 r->iv[i + r->in].op = op;
148 r->iv[i + r->in].argx = w;
149 r->iv[i + r->in].argy = b;
150 }
151 }
152 DA_DESTROY(&iv);
153 }
154
155 /* --- @mpreduce_destroy@ --- *
156 *
157 * Arguments: @mpreduce *r@ = structure to free
158 *
159 * Returns: ---
160 *
161 * Use: Reclaims the resources from a reduction context.
162 */
163
164 void mpreduce_destroy(mpreduce *r)
165 {
166 mp_drop(r->p);
167 xfree(r->iv);
168 }
169
170 /* --- @mpreduce_dump@ --- *
171 *
172 * Arguments: @mpreduce *r@ = structure to dump
173 * @FILE *fp@ = file to dump on
174 *
175 * Returns: ---
176 *
177 * Use: Dumps a reduction context.
178 */
179
180 void mpreduce_dump(mpreduce *r, FILE *fp)
181 {
182 size_t i;
183 static const char *opname[] = { "add", "addshift", "sub", "subshift" };
184
185 fprintf(fp, "mod = "); mp_writefile(r->p, fp, 16);
186 fprintf(fp, "\n lim = %lu; s = %d\n", (unsigned long)r->lim, r->s);
187 for (i = 0; i < r->in; i++) {
188 assert(r->iv[i].op < N(opname));
189 fprintf(fp, " %s %lu %lu\n",
190 opname[r->iv[i].op],
191 (unsigned long)r->iv[i].argx,
192 (unsigned long)r->iv[i].argy);
193 }
194 if (r->s) {
195 fprintf(fp, "tail end charlie\n");
196 for (i = r->in; i < 2 * r->in; i++) {
197 assert(r->iv[i].op < N(opname));
198 fprintf(fp, " %s %lu %lu\n",
199 opname[r->iv[i].op],
200 (unsigned long)r->iv[i].argx,
201 (unsigned long)r->iv[i].argy);
202 }
203 }
204 }
205
206 /* --- @mpreduce_do@ --- *
207 *
208 * Arguments: @mpreduce *r@ = reduction context
209 * @mp *d@ = destination
210 * @mp *x@ = source
211 *
212 * Returns: Destination, @x@ reduced modulo the reduction poly.
213 */
214
215 static void run(const mpreduce_instr *i, const mpreduce_instr *il,
216 mpw *v, mpw z)
217 {
218 for (; i < il; i++) {
219 #ifdef DEBUG
220 mp vv;
221 mp_build(&vv, v - i->argx, v + 1);
222 printf(" 0x"); mp_writefile(&vv, stdout, 16);
223 printf(" %c (0x%lx << %u) == 0x",
224 (i->op & ~1u) == MPRI_ADD ? '+' : '-',
225 (unsigned long)z,
226 i->argy);
227 #endif
228 switch (i->op) {
229 case MPRI_ADD: MPX_UADDN(v - i->argx, v + 1, z); break;
230 case MPRI_ADDLSL: mpx_uaddnlsl(v - i->argx, v + 1, z, i->argy); break;
231 case MPRI_SUB: MPX_USUBN(v - i->argx, v + 1, z); break;
232 case MPRI_SUBLSL: mpx_usubnlsl(v - i->argx, v + 1, z, i->argy); break;
233 default:
234 abort();
235 }
236 #ifdef DEBUG
237 mp_build(&vv, v - i->argx, v + 1);
238 mp_writefile(&vv, stdout, 16);
239 printf("\n");
240 #endif
241 }
242 }
243
244 mp *mpreduce_do(mpreduce *r, mp *d, mp *x)
245 {
246 mpw *v, *vl;
247 const mpreduce_instr *il;
248 mpw z;
249
250 #ifdef DEBUG
251 mp *_r = 0, *_rr = 0;
252 #endif
253
254 /* --- If source is negative, divide --- */
255
256 if (MP_ISNEG(x)) {
257 mp_div(0, &d, x, r->p);
258 return (d);
259 }
260
261 /* --- Try to reuse the source's space --- */
262
263 MP_COPY(x);
264 if (d) MP_DROP(d);
265 MP_DEST(x, MP_LEN(x), x->f);
266
267 /* --- Do the reduction --- */
268
269 #ifdef DEBUG
270 _r = MP_NEW;
271 mp_div(0, &_r, x, r->p);
272 MP_PRINTX("x", x);
273 _rr = 0;
274 #endif
275
276 il = r->iv + r->in;
277 if (MP_LEN(x) >= r->lim) {
278 v = x->v + r->lim;
279 vl = x->vl;
280 while (vl-- > v) {
281 while (*vl) {
282 z = *vl;
283 *vl = 0;
284 run(r->iv, il, vl, z);
285 #ifdef DEBUG
286 MP_PRINTX("x", x);
287 mp_div(0, &_rr, x, r->p);
288 assert(MP_EQ(_r, _rr));
289 #endif
290 }
291 }
292 if (r->s) {
293 while (*vl >> r->s) {
294 z = *vl >> r->s;
295 *vl &= ((1 << r->s) - 1);
296 run(r->iv + r->in, il + r->in, vl, z);
297 #ifdef DEBUG
298 MP_PRINTX("x", x);
299 mp_div(0, &_rr, x, r->p);
300 assert(MP_EQ(_r, _rr));
301 #endif
302 }
303 }
304 }
305
306 /* --- Finishing touches --- */
307
308 MP_SHRINK(x);
309 if (MP_CMP(x, >=, r->p))
310 x = mp_sub(x, x, r->p);
311
312 /* --- Done --- */
313
314 #ifdef DEBUG
315 assert(MP_EQ(_r, x));
316 mp_drop(_r);
317 mp_drop(_rr);
318 #endif
319 return (x);
320 }
321
322 /* --- @mpreduce_exp@ --- *
323 *
324 * Arguments: @mpreduce *mr@ = pointer to reduction context
325 * @mp *d@ = fake destination
326 * @mp *a@ = base
327 * @mp *e@ = exponent
328 *
329 * Returns: Result, %$a^e \bmod m$%.
330 */
331
332 mp *mpreduce_exp(mpreduce *mr, mp *d, mp *a, mp *e)
333 {
334 mp *x = MP_ONE;
335 mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW;
336
337 MP_SHRINK(e);
338 if (!MP_LEN(e))
339 ;
340 else if (MP_LEN(e) < EXP_THRESH)
341 EXP_SIMPLE(x, a, e);
342 else
343 EXP_WINDOW(x, a, e);
344 mp_drop(d);
345 mp_drop(spare);
346 return (x);
347 }
348
349 /*----- Test rig ----------------------------------------------------------*/
350
351
352 #ifdef TEST_RIG
353
354 #define MP(x) mp_readstring(MP_NEW, #x, 0, 0)
355
356 static int vreduce(dstr *v)
357 {
358 mp *d = *(mp **)v[0].buf;
359 mp *n = *(mp **)v[1].buf;
360 mp *r = *(mp **)v[2].buf;
361 mp *c;
362 int ok = 1;
363 mpreduce rr;
364
365 mpreduce_create(&rr, d);
366 c = mpreduce_do(&rr, MP_NEW, n);
367 if (!MP_EQ(c, r)) {
368 fprintf(stderr, "\n*** reduction failed\n*** ");
369 mpreduce_dump(&rr, stderr);
370 fprintf(stderr, "\n*** n = "); mp_writefile(n, stderr, 10);
371 fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 10);
372 fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 10);
373 fprintf(stderr, "\n");
374 ok = 0;
375 }
376 mpreduce_destroy(&rr);
377 mp_drop(n); mp_drop(d); mp_drop(r); mp_drop(c);
378 assert(mparena_count(MPARENA_GLOBAL) == 0);
379 return (ok);
380 }
381
382 static int vmodexp(dstr *v)
383 {
384 mp *p = *(mp **)v[0].buf;
385 mp *g = *(mp **)v[1].buf;
386 mp *x = *(mp **)v[2].buf;
387 mp *r = *(mp **)v[3].buf;
388 mp *c;
389 int ok = 1;
390 mpreduce rr;
391
392 mpreduce_create(&rr, p);
393 c = mpreduce_exp(&rr, MP_NEW, g, x);
394 if (!MP_EQ(c, r)) {
395 fprintf(stderr, "\n*** modexp failed\n*** ");
396 fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 10);
397 fprintf(stderr, "\n*** g = "); mp_writefile(g, stderr, 10);
398 fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 10);
399 fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 10);
400 fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 10);
401 fprintf(stderr, "\n");
402 ok = 0;
403 }
404 mpreduce_destroy(&rr);
405 mp_drop(p); mp_drop(g); mp_drop(r); mp_drop(x); mp_drop(c);
406 assert(mparena_count(MPARENA_GLOBAL) == 0);
407 return (ok);
408 }
409
410 static test_chunk defs[] = {
411 { "reduce", vreduce, { &type_mp, &type_mp, &type_mp, 0 } },
412 { "modexp", vmodexp, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
413 { 0, 0, { 0 } }
414 };
415
416 int main(int argc, char *argv[])
417 {
418 test_run(argc, argv, defs, SRCDIR"/tests/mpreduce");
419 return (0);
420 }
421
422 #endif
423
424 /*----- That's all, folks -------------------------------------------------*/
Catacomb cryptographic library (downstream of http://git.distorted.org.uk/~mdw/catacomb/)