projects / u / mdw / catacomb / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
Change interface for suggested destinations.
[u/mdw/catacomb] / mpmont.c
1 /* -*-c-*-
2 *
3 * $Id: mpmont.c,v 1.6 1999/12/10 23:18:39 mdw Exp $
4 *
5 * Montgomery reduction
6 *
7 * (c) 1999 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: mpmont.c,v $
33 * Revision 1.6 1999/12/10 23:18:39 mdw
34 * Change interface for suggested destinations.
35 *
36 * Revision 1.5 1999/11/22 13:58:40 mdw
37 * Add an option to disable Montgomery reduction, so that performance
38 * comparisons can be done.
39 *
40 * Revision 1.4 1999/11/21 12:27:06 mdw
41 * Remove a division from the Montgomery setup by calculating
42 * %$R^2 \bmod m$% first and then %$R \bmod m$% by Montgomery reduction of
43 * %$R^2$%.
44 *
45 * Revision 1.3 1999/11/21 11:35:10 mdw
46 * Performance improvement: use @mp_sqr@ and @mpmont_reduce@ instead of
47 * @mpmont_mul@ for squaring in exponentiation.
48 *
49 * Revision 1.2 1999/11/19 13:17:26 mdw
50 * Add extra interface to exponentiation which returns a Montgomerized
51 * result.
52 *
53 * Revision 1.1 1999/11/17 18:02:16 mdw
54 * New multiprecision integer arithmetic suite.
55 *
56 */
57
58 /*----- Header files ------------------------------------------------------*/
59
60 #include "mp.h"
61 #include "mpmont.h"
62
63 /*----- Tweakables --------------------------------------------------------*/
64
65 /* --- @MPMONT_DISABLE@ --- *
66 *
67 * Replace all the clever Montgomery reduction with good old-fashioned long
68 * division.
69 */
70
71 /* #define MPMONT_DISABLE */
72
73 /*----- Main code ---------------------------------------------------------*/
74
75 /* --- @mpmont_create@ --- *
76 *
77 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
78 * @mp *m@ = modulus to use
79 *
80 * Returns: ---
81 *
82 * Use: Initializes a Montgomery reduction context ready for use.
83 * The argument @m@ must be a positive odd integer.
84 */
85
86 #ifdef MPMONT_DISABLE
87
88 void mpmont_create(mpmont *mm, mp *m)
89 {
90 mp_shrink(m);
91 mm->m = MP_COPY(m);
92 mm->r = MP_ONE;
93 mm->r2 = MP_ONE;
94 }
95
96 #else
97
98 void mpmont_create(mpmont *mm, mp *m)
99 {
100 /* --- Validate the arguments --- */
101
102 assert(((void)"Montgomery modulus must be positive",
103 (m->f & MP_NEG) == 0));
104 assert(((void)"Montgomery modulus must be odd", m->v[0] & 1));
105
106 /* --- Take a copy of the modulus --- */
107
108 mp_shrink(m);
109 mm->m = MP_COPY(m);
110
111 /* --- Find the magic value @mi@ --- *
112 *
113 * This is a slightly grungy way of solving the problem, but it does work.
114 */
115
116 {
117 mpw av[2] = { 0, 1 };
118 mp a, b;
119 mp *i = MP_NEW;
120 mpw mi;
121
122 mp_build(&a, av, av + 2);
123 mp_build(&b, m->v, m->v + 1);
124 mp_gcd(0, 0, &i, &a, &b);
125 mi = i->v[0];
126 if (!(i->f & MP_NEG))
127 mi = MPW(-mi);
128 mm->mi = mi;
129 MP_DROP(i);
130 }
131
132 /* --- Discover the values %$R \bmod m$% and %$R^2 \bmod m$% --- */
133
134 {
135 size_t l = MP_LEN(m);
136 mp *r = mp_create(2 * l + 1);
137
138 mm->shift = l * MPW_BITS;
139 MPX_ZERO(r->v, r->vl - 1);
140 r->vl[-1] = 1;
141 mm->r2 = MP_NEW;
142 mp_div(0, &mm->r2, r, m);
143 mm->r = mpmont_reduce(mm, MP_NEW, mm->r2);
144 MP_DROP(r);
145 }
146 }
147
148 #endif
149
150 /* --- @mpmont_destroy@ --- *
151 *
152 * Arguments: @mpmont *mm@ = pointer to a Montgomery reduction context
153 *
154 * Returns: ---
155 *
156 * Use: Disposes of a context when it's no longer of any use to
157 * anyone.
158 */
159
160 void mpmont_destroy(mpmont *mm)
161 {
162 MP_DROP(mm->m);
163 MP_DROP(mm->r);
164 MP_DROP(mm->r2);
165 }
166
167 /* --- @mpmont_reduce@ --- *
168 *
169 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
170 * @mp *d@ = destination
171 * @mp *a@ = source, assumed positive
172 *
173 * Returns: Result, %$a R^{-1} \bmod m$%.
174 */
175
176 #ifdef MPMONT_DISABLE
177
178 mp *mpmont_reduce(mpmont *mm, mp *d, mp *a)
179 {
180 mp_div(0, &d, a, mm->m);
181 return (d);
182 }
183
184 #else
185
186 mp *mpmont_reduce(mpmont *mm, mp *d, mp *a)
187 {
188 mpw *dv, *dvl;
189 mpw *mv, *mvl;
190 size_t n;
191
192 /* --- Initial conditioning of the arguments --- */
193
194 n = MP_LEN(mm->m);
195
196 if (d == a)
197 MP_MODIFY(d, 2 * n + 1);
198 else {
199 MP_MODIFY(d, 2 * n + 1);
200 memcpy(d->v, a->v, MPWS(MP_LEN(a)));
201 memset(d->v + MP_LEN(a), 0, MPWS(MP_LEN(d) - MP_LEN(a)));
202 }
203
204 dv = d->v; dvl = d->vl;
205 mv = mm->m->v; mvl = mm->m->vl;
206
207 /* --- Let's go to work --- */
208
209 while (n--) {
210 mpw u = MPW(*dv * mm->mi);
211 MPX_UMLAN(dv, dvl, mv, mvl, u);
212 dv++;
213 }
214
215 /* --- Done --- */
216
217 memmove(d->v, dv, MPWS(dvl - dv));
218 d->vl -= dv - d->v;
219 MP_SHRINK(d);
220 d->f = a->f & MP_BURN;
221 if (MP_CMP(d, >=, mm->m))
222 d = mp_sub(d, d, mm->m);
223 return (d);
224 }
225
226 #endif
227
228 /* --- @mpmont_mul@ --- *
229 *
230 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
231 * @mp *d@ = destination
232 * @mp *a, *b@ = sources, assumed positive
233 *
234 * Returns: Result, %$a b R^{-1} \bmod m$%.
235 */
236
237 #ifdef MPMONT_DISABLE
238
239 mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b)
240 {
241 d = mp_mul(d, a, b);
242 mp_div(0, &d, d, mm->m);
243 return (d);
244 }
245
246 #else
247
248 mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b)
249 {
250 if (MP_LEN(a) > KARATSUBA_CUTOFF && MP_LEN(b) > KARATSUBA_CUTOFF) {
251 d = mp_mul(d, a, b);
252 d = mpmont_reduce(mm, d, d);
253 } else {
254 mpw *dv, *dvl;
255 mpw *av, *avl;
256 mpw *bv, *bvl;
257 mpw *mv, *mvl;
258 mpw y;
259 size_t n, i;
260
261 /* --- Initial conditioning of the arguments --- */
262
263 if (MP_LEN(a) > MP_LEN(b)) {
264 mp *t = a; a = b; b = t;
265 }
266 n = MP_LEN(mm->m);
267
268 a = MP_COPY(a);
269 b = MP_COPY(b);
270 MP_MODIFY(d, 2 * n + 1);
271 dv = d->v; dvl = d->vl;
272 MPX_ZERO(dv, dvl);
273 av = a->v; avl = a->vl;
274 bv = b->v; bvl = b->vl;
275 mv = mm->m->v; mvl = mm->m->vl;
276 y = *bv;
277
278 /* --- Montgomery multiplication phase --- */
279
280 i = 0;
281 while (i < n && av < avl) {
282 mpw x = *av++;
283 mpw u = MPW((*dv + x * y) * mm->mi);
284 MPX_UMLAN(dv, dvl, bv, bvl, x);
285 MPX_UMLAN(dv, dvl, mv, mvl, u);
286 dv++;
287 i++;
288 }
289
290 /* --- Simpler Montgomery reduction phase --- */
291
292 while (i < n) {
293 mpw u = MPW(*dv * mm->mi);
294 MPX_UMLAN(dv, dvl, mv, mvl, u);
295 dv++;
296 i++;
297 }
298
299 /* --- Done --- */
300
301 memmove(d->v, dv, MPWS(dvl - dv));
302 d->vl -= dv - d->v;
303 MP_SHRINK(d);
304 d->f = (a->f | b->f) & MP_BURN;
305 if (MP_CMP(d, >=, mm->m))
306 d = mp_sub(d, d, mm->m);
307 MP_DROP(a);
308 MP_DROP(b);
309 }
310
311 return (d);
312 }
313
314 #endif
315
316 /* --- @mpmont_expr@ --- *
317 *
318 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
319 * @mp *d@ = fake destination
320 * @mp *a@ = base
321 * @mp *e@ = exponent
322 *
323 * Returns: Result, %$a^e R \bmod m$%.
324 */
325
326 mp *mpmont_expr(mpmont *mm, mp *d, mp *a, mp *e)
327 {
328 mpscan sc;
329 mp *ar = mpmont_mul(mm, MP_NEW, a, mm->r2);
330 mp *x = MP_COPY(mm->r);
331 mp *spare = MP_NEW;
332
333 mp_scan(&sc, e);
334
335 if (MP_STEP(&sc)) {
336 size_t sq = 0;
337 for (;;) {
338 mp *dd;
339 if (MP_BIT(&sc)) {
340 while (sq) {
341 dd = mp_sqr(spare, ar);
342 dd = mpmont_reduce(mm, dd, dd);
343 spare = ar; ar = dd;
344 sq--;
345 }
346 dd = mpmont_mul(mm, spare, x, ar);
347 spare = x; x = dd;
348 }
349 sq++;
350 if (!MP_STEP(&sc))
351 break;
352 }
353 }
354 MP_DROP(ar);
355 if (spare != MP_NEW)
356 MP_DROP(spare);
357 if (d != MP_NEW)
358 MP_DROP(d);
359 return (x);
360 }
361
362 /* --- @mpmont_exp@ --- *
363 *
364 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
365 * @mp *d@ = fake destination
366 * @mp *a@ = base
367 * @mp *e@ = exponent
368 *
369 * Returns: Result, %$a^e \bmod m$%.
370 */
371
372 mp *mpmont_exp(mpmont *mm, mp *d, mp *a, mp *e)
373 {
374 d = mpmont_expr(mm, d, a, e);
375 d = mpmont_reduce(mm, d, d);
376 return (d);
377 }
378
379 /*----- Test rig ----------------------------------------------------------*/
380
381 #ifdef TEST_RIG
382
383 static int tcreate(dstr *v)
384 {
385 mp *m = *(mp **)v[0].buf;
386 mp *mi = *(mp **)v[1].buf;
387 mp *r = *(mp **)v[2].buf;
388 mp *r2 = *(mp **)v[3].buf;
389
390 mpmont mm;
391 int ok = 1;
392
393 mpmont_create(&mm, m);
394
395 if (mm.mi != mi->v[0]) {
396 fprintf(stderr, "\n*** bad mi: found %lu, expected %lu",
397 (unsigned long)mm.mi, (unsigned long)mi->v[0]);
398 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
399 fputc('\n', stderr);
400 ok = 0;
401 }
402
403 if (MP_CMP(mm.r, !=, r)) {
404 fputs("\n*** bad r", stderr);
405 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
406 fputs("\nexpected ", stderr); mp_writefile(r, stderr, 10);
407 fputs("\n found ", stderr); mp_writefile(mm.r, stderr, 10);
408 fputc('\n', stderr);
409 ok = 0;
410 }
411
412 if (MP_CMP(mm.r2, !=, r2)) {
413 fputs("\n*** bad r2", stderr);
414 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
415 fputs("\nexpected ", stderr); mp_writefile(r2, stderr, 10);
416 fputs("\n found ", stderr); mp_writefile(mm.r2, stderr, 10);
417 fputc('\n', stderr);
418 ok = 0;
419 }
420
421 MP_DROP(m);
422 MP_DROP(mi);
423 MP_DROP(r);
424 MP_DROP(r2);
425 mpmont_destroy(&mm);
426 assert(mparena_count(MPARENA_GLOBAL) == 0);
427 return (ok);
428 }
429
430 static int tmul(dstr *v)
431 {
432 mp *m = *(mp **)v[0].buf;
433 mp *a = *(mp **)v[1].buf;
434 mp *b = *(mp **)v[2].buf;
435 mp *r = *(mp **)v[3].buf;
436 int ok = 1;
437
438 mpmont mm;
439 mpmont_create(&mm, m);
440
441 {
442 mp *qr = mp_mul(MP_NEW, a, b);
443 mp_div(0, &qr, qr, m);
444
445 if (MP_CMP(qr, !=, r)) {
446 fputs("\n*** classical modmul failed", stderr);
447 fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
448 fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
449 fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
450 fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
451 fputs("\nqr = ", stderr); mp_writefile(qr, stderr, 10);
452 fputc('\n', stderr);
453 ok = 0;
454 }
455
456 mp_drop(qr);
457 }
458
459 {
460 mp *ar = mpmont_mul(&mm, MP_NEW, a, mm.r2);
461 mp *br = mpmont_mul(&mm, MP_NEW, b, mm.r2);
462 mp *mr = mpmont_mul(&mm, MP_NEW, ar, br);
463 mr = mpmont_reduce(&mm, mr, mr);
464 if (MP_CMP(mr, !=, r)) {
465 fputs("\n*** montgomery modmul failed", stderr);
466 fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
467 fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
468 fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
469 fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
470 fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10);
471 fputc('\n', stderr);
472 ok = 0;
473 }
474 MP_DROP(ar); MP_DROP(br);
475 mp_drop(mr);
476 }
477
478
479 MP_DROP(m);
480 MP_DROP(a);
481 MP_DROP(b);
482 MP_DROP(r);
483 mpmont_destroy(&mm);
484 assert(mparena_count(MPARENA_GLOBAL) == 0);
485 return ok;
486 }
487
488 static int texp(dstr *v)
489 {
490 mp *m = *(mp **)v[0].buf;
491 mp *a = *(mp **)v[1].buf;
492 mp *b = *(mp **)v[2].buf;
493 mp *r = *(mp **)v[3].buf;
494 mp *mr;
495 int ok = 1;
496
497 mpmont mm;
498 mpmont_create(&mm, m);
499
500 mr = mpmont_exp(&mm, MP_NEW, a, b);
501
502 if (MP_CMP(mr, !=, r)) {
503 fputs("\n*** montgomery modexp failed", stderr);
504 fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
505 fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
506 fputs("\n e = ", stderr); mp_writefile(b, stderr, 10);
507 fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
508 fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10);
509 fputc('\n', stderr);
510 ok = 0;
511 }
512
513 MP_DROP(m);
514 MP_DROP(a);
515 MP_DROP(b);
516 MP_DROP(r);
517 MP_DROP(mr);
518 mpmont_destroy(&mm);
519 assert(mparena_count(MPARENA_GLOBAL) == 0);
520 return ok;
521 }
522
523
524 static test_chunk tests[] = {
525 { "create", tcreate, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
526 { "mul", tmul, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
527 { "exp", texp, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
528 { 0, 0, { 0 } },
529 };
530
531 int main(int argc, char *argv[])
532 {
533 sub_init();
534 test_run(argc, argv, tests, SRCDIR "/tests/mpmont");
535 return (0);
536 }
537
538 #endif
539
540 /*----- That's all, folks -------------------------------------------------*/
Catacomb cryptographic library (downstream of http://git.distorted.org.uk/~mdw/catacomb/)