projects / u / mdw / catacomb / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Change interface for suggested destinations.
[u/mdw/catacomb] / mpmont.c
CommitLineData
d3409d5e 1/* -*-c-*-
2 *
ef5f4810 3 * $Id: mpmont.c,v 1.6 1999/12/10 23:18:39 mdw Exp $
d3409d5e 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 $
ef5f4810 33 * Revision 1.6 1999/12/10 23:18:39 mdw
34 * Change interface for suggested destinations.
35 *
52e4b041 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 *
93feaa6e 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 *
79a34029 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 *
17ad212e 49 * Revision 1.2 1999/11/19 13:17:26 mdw
50 * Add extra interface to exponentiation which returns a Montgomerized
51 * result.
52 *
d3409d5e 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
52e4b041 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
d3409d5e 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.
ef5f4810 83 * The argument @m@ must be a positive odd integer.
d3409d5e 84 */
85
52e4b041 86#ifdef MPMONT_DISABLE
87
88void 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
d3409d5e 98void mpmont_create(mpmont *mm, mp *m)
99{
ef5f4810 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
d3409d5e 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;
ef5f4810 119 mp *i = MP_NEW;
d3409d5e 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);
93feaa6e 136 mp *r = mp_create(2 * l + 1);
d3409d5e 137
138 mm->shift = l * MPW_BITS;
139 MPX_ZERO(r->v, r->vl - 1);
140 r->vl[-1] = 1;
93feaa6e 141 mm->r2 = MP_NEW;
d3409d5e 142 mp_div(0, &mm->r2, r, m);
93feaa6e 143 mm->r = mpmont_reduce(mm, MP_NEW, mm->r2);
d3409d5e 144 MP_DROP(r);
145 }
146}
147
52e4b041 148#endif
149
d3409d5e 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
160void 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
ef5f4810 171 * @mp *a@ = source, assumed positive
d3409d5e 172 *
173 * Returns: Result, %$a R^{-1} \bmod m$%.
174 */
175
52e4b041 176#ifdef MPMONT_DISABLE
177
ef5f4810 178mp *mpmont_reduce(mpmont *mm, mp *d, mp *a)
52e4b041 179{
180 mp_div(0, &d, a, mm->m);
181 return (d);
182}
183
184#else
185
ef5f4810 186mp *mpmont_reduce(mpmont *mm, mp *d, mp *a)
d3409d5e 187{
188 mpw *dv, *dvl;
ef5f4810 189 mpw *mv, *mvl;
d3409d5e 190 size_t n;
191
192 /* --- Initial conditioning of the arguments --- */
193
194 n = MP_LEN(mm->m);
195
196 if (d == a)
79a34029 197 MP_MODIFY(d, 2 * n + 1);
d3409d5e 198 else {
79a34029 199 MP_MODIFY(d, 2 * n + 1);
d3409d5e 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;
17ad212e 221 if (MP_CMP(d, >=, mm->m))
222 d = mp_sub(d, d, mm->m);
d3409d5e 223 return (d);
224}
225
52e4b041 226#endif
227
d3409d5e 228/* --- @mpmont_mul@ --- *
229 *
230 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
231 * @mp *d@ = destination
ef5f4810 232 * @mp *a, *b@ = sources, assumed positive
d3409d5e 233 *
234 * Returns: Result, %$a b R^{-1} \bmod m$%.
235 */
236
52e4b041 237#ifdef MPMONT_DISABLE
238
ef5f4810 239mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b)
52e4b041 240{
241 d = mp_mul(d, a, b);
242 mp_div(0, &d, d, mm->m);
243 return (d);
244}
245
246#else
247
ef5f4810 248mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b)
d3409d5e 249{
ef5f4810 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);
d3409d5e 267
ef5f4810 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 }
d3409d5e 289
ef5f4810 290 /* --- Simpler Montgomery reduction phase --- */
d3409d5e 291
ef5f4810 292 while (i < n) {
293 mpw u = MPW(*dv * mm->mi);
294 MPX_UMLAN(dv, dvl, mv, mvl, u);
295 dv++;
296 i++;
297 }
d3409d5e 298
ef5f4810 299 /* --- Done --- */
d3409d5e 300
ef5f4810 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);
d3409d5e 309 }
310
d3409d5e 311 return (d);
312}
313
52e4b041 314#endif
315
17ad212e 316/* --- @mpmont_expr@ --- *
d3409d5e 317 *
318 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
ef5f4810 319 * @mp *d@ = fake destination
320 * @mp *a@ = base
321 * @mp *e@ = exponent
d3409d5e 322 *
17ad212e 323 * Returns: Result, %$a^e R \bmod m$%.
d3409d5e 324 */
325
ef5f4810 326mp *mpmont_expr(mpmont *mm, mp *d, mp *a, mp *e)
d3409d5e 327{
328 mpscan sc;
329 mp *ar = mpmont_mul(mm, MP_NEW, a, mm->r2);
ef5f4810 330 mp *x = MP_COPY(mm->r);
17ad212e 331 mp *spare = MP_NEW;
d3409d5e 332
333 mp_scan(&sc, e);
334
335 if (MP_STEP(&sc)) {
17ad212e 336 size_t sq = 0;
d3409d5e 337 for (;;) {
338 mp *dd;
339 if (MP_BIT(&sc)) {
17ad212e 340 while (sq) {
79a34029 341 dd = mp_sqr(spare, ar);
342 dd = mpmont_reduce(mm, dd, dd);
17ad212e 343 spare = ar; ar = dd;
344 sq--;
345 }
ef5f4810 346 dd = mpmont_mul(mm, spare, x, ar);
347 spare = x; x = dd;
d3409d5e 348 }
17ad212e 349 sq++;
d3409d5e 350 if (!MP_STEP(&sc))
351 break;
d3409d5e 352 }
353 }
354 MP_DROP(ar);
17ad212e 355 if (spare != MP_NEW)
356 MP_DROP(spare);
ef5f4810 357 if (d != MP_NEW)
358 MP_DROP(d);
359 return (x);
17ad212e 360}
361
362/* --- @mpmont_exp@ --- *
363 *
364 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
ef5f4810 365 * @mp *d@ = fake destination
366 * @mp *a@ = base
367 * @mp *e@ = exponent
17ad212e 368 *
369 * Returns: Result, %$a^e \bmod m$%.
370 */
371
ef5f4810 372mp *mpmont_exp(mpmont *mm, mp *d, mp *a, mp *e)
17ad212e 373{
ef5f4810 374 d = mpmont_expr(mm, d, a, e);
17ad212e 375 d = mpmont_reduce(mm, d, d);
376 return (d);
d3409d5e 377}
378
379/*----- Test rig ----------------------------------------------------------*/
380
381#ifdef TEST_RIG
382
383static 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);
17ad212e 407 fputs("\n found ", stderr); mp_writefile(mm.r, stderr, 10);
d3409d5e 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);
17ad212e 416 fputs("\n found ", stderr); mp_writefile(mm.r2, stderr, 10);
d3409d5e 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);
ef5f4810 426 assert(mparena_count(MPARENA_GLOBAL) == 0);
d3409d5e 427 return (ok);
428}
429
430static 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;
d3409d5e 436 int ok = 1;
437
438 mpmont mm;
439 mpmont_create(&mm, m);
440
441 {
79a34029 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 {
d3409d5e 460 mp *ar = mpmont_mul(&mm, MP_NEW, a, mm.r2);
461 mp *br = mpmont_mul(&mm, MP_NEW, b, mm.r2);
79a34029 462 mp *mr = mpmont_mul(&mm, MP_NEW, ar, br);
d3409d5e 463 mr = mpmont_reduce(&mm, mr, mr);
79a34029 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 }
d3409d5e 474 MP_DROP(ar); MP_DROP(br);
79a34029 475 mp_drop(mr);
d3409d5e 476 }
477
d3409d5e 478
479 MP_DROP(m);
480 MP_DROP(a);
481 MP_DROP(b);
482 MP_DROP(r);
d3409d5e 483 mpmont_destroy(&mm);
ef5f4810 484 assert(mparena_count(MPARENA_GLOBAL) == 0);
d3409d5e 485 return ok;
486}
487
488static 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
ef5f4810 500 mr = mpmont_exp(&mm, MP_NEW, a, b);
d3409d5e 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);
ef5f4810 519 assert(mparena_count(MPARENA_GLOBAL) == 0);
d3409d5e 520 return ok;
521}
522
523
524static test_chunk tests[] = {
ef5f4810 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 } },
d3409d5e 528 { 0, 0, { 0 } },
529};
530
531int 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/)