base/asm-common.h: Add some debugging macros.
[catacomb] / math / mpmont.c
CommitLineData
d3409d5e 1/* -*-c-*-
2 *
d3409d5e 3 * Montgomery reduction
4 *
5 * (c) 1999 Straylight/Edgeware
6 */
7
45c0fd36 8/*----- Licensing notice --------------------------------------------------*
d3409d5e 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.
45c0fd36 16 *
d3409d5e 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.
45c0fd36 21 *
d3409d5e 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
d3409d5e 28/*----- Header files ------------------------------------------------------*/
29
30#include "mp.h"
31#include "mpmont.h"
32
52e4b041 33/*----- Tweakables --------------------------------------------------------*/
34
35/* --- @MPMONT_DISABLE@ --- *
36 *
37 * Replace all the clever Montgomery reduction with good old-fashioned long
38 * division.
39 */
40
41/* #define MPMONT_DISABLE */
42
0e70bd46
MW
43/*----- Low-level implementation ------------------------------------------*/
44
45#ifndef MPMONT_DISABLE
46
47/* --- @redccore@ --- *
48 *
49 * Arguments: @mpw *dv, *dvl@ = base and limit of source/destination
50 * @const mpw *mv@ = base of modulus %$m$%
51 * @size_t n@ = length of modulus
52 * @const mpw *mi@ = base of REDC coefficient %$m'$%
53 *
54 * Returns: ---
55 *
56 * Use: Let %$a$% be the input operand. Store in %$d$% the value
57 * %$a + (m' a \bmod R) m$%. The destination has space for at
58 * least %$2 n + 1$% words of result.
59 */
60
61static void redccore(mpw *dv, mpw *dvl, const mpw *mv,
62 size_t n, const mpw *mi)
63{
64 mpw mi0 = *mi;
65 size_t i;
66
67 for (i = 0; i < n; i++) {
68 MPX_UMLAN(dv, dvl, mv, mv + n, MPW(*dv*mi0));
69 dv++;
70 }
71}
72
73/* --- @redccore@ --- *
74 *
75 * Arguments: @mpw *dv, *dvl@ = base and limit of source/destination
76 * @const mpw *av, *avl@ = base and limit of first multiplicand
77 * @const mpw *bv, *bvl@ = base and limit of second multiplicand
78 * @const mpw *mv@ = base of modulus %$m$%
79 * @size_t n@ = length of modulus
80 * @const mpw *mi@ = base of REDC coefficient %$m'$%
81 *
82 * Returns: ---
83 *
84 * Use: Let %$a$% and %$b$% be the multiplicands. Let %$w = a b$%.
85 * Store in %$d$% the value %$a b + (m' a b \bmod R) m$%.
86 */
87
88static void mulcore(mpw *dv, mpw *dvl,
89 const mpw *av, const mpw *avl,
90 const mpw *bv, const mpw *bvl,
91 const mpw *mv, size_t n, const mpw *mi)
92{
93 mpw ai, b0, y, mi0 = *mi;
94 const mpw *tv, *tvl;
95 const mpw *mvl = mv + n;
96 size_t i = 0;
97
98 /* --- Initial setup --- */
99
100 MPX_ZERO(dv, dvl);
101 if (avl - av > bvl - bv) {
102 tv = av; av = bv; bv = tv;
103 tvl = avl; avl = bvl; bvl = tvl;
104 }
105 b0 = *bv;
106
107 /* --- Multiply, until we run out of multiplicand --- */
108
109 while (i < n && av < avl) {
110 ai = *av++;
111 y = MPW((*dv + ai*b0)*mi0);
112 MPX_UMLAN(dv, dvl, bv, bvl, ai);
113 MPX_UMLAN(dv, dvl, mv, mvl, y);
114 dv++; i++;
115 }
116
117 /* --- Continue reducing until we run out of modulus --- */
118
119 while (i < n) {
120 y = MPW(*dv*mi0);
121 MPX_UMLAN(dv, dvl, mv, mvl, y);
122 dv++; i++;
123 }
124}
125
126/* --- @finish@ --- *
127 *
128 * Arguments: @mpmont *mm@ = pointer to a Montgomery reduction context
129 * *mp *d@ = pointer to mostly-reduced operand
130 *
131 * Returns: ---
132 *
133 * Use: Applies the finishing touches to Montgomery reduction. The
134 * operand @d@ is a multiple of %$R%$ at this point, so it needs
135 * to be shifted down; the result might need a further
136 * subtraction to get it into the right interval; and we may
137 * need to do an additional subtraction if %$d$% is negative.
138 */
139
140static void finish(mpmont *mm, mp *d)
141{
142 mpw *dv = d->v, *dvl = d->vl;
143 size_t n = mm->n;
144
145 memmove(dv, dv + n, MPWS(dvl - (dv + n)));
146 dvl -= n;
147
148 if (MPX_UCMP(dv, dvl, >=, mm->m->v, mm->m->vl))
149 mpx_usub(dv, dvl, dv, dvl, mm->m->v, mm->m->vl);
150
151 if (d->f & MP_NEG) {
152 mpx_usub(dv, dvl, mm->m->v, mm->m->vl, dv, dvl);
153 d->f &= ~MP_NEG;
154 }
155
156 d->vl = dvl;
157 MP_SHRINK(d);
158}
159
160#endif
161
4640a0dd 162/*----- Reduction and multiplication --------------------------------------*/
d3409d5e 163
164/* --- @mpmont_create@ --- *
165 *
166 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
167 * @mp *m@ = modulus to use
168 *
f4535c64 169 * Returns: Zero on success, nonzero on error.
d3409d5e 170 *
171 * Use: Initializes a Montgomery reduction context ready for use.
ef5f4810 172 * The argument @m@ must be a positive odd integer.
d3409d5e 173 */
174
52e4b041 175#ifdef MPMONT_DISABLE
176
f4535c64 177int mpmont_create(mpmont *mm, mp *m)
52e4b041 178{
179 mp_shrink(m);
180 mm->m = MP_COPY(m);
181 mm->r = MP_ONE;
182 mm->r2 = MP_ONE;
f5f35081 183 mm->mi = MP_ONE;
f4535c64 184 return (0);
52e4b041 185}
186
187#else
188
f4535c64 189int mpmont_create(mpmont *mm, mp *m)
d3409d5e 190{
f5f35081 191 size_t n = MP_LEN(m);
d34decd2 192 mp *r2 = mp_new(2 * n + 1, 0);
f5f35081 193 mp r;
194
d3409d5e 195 /* --- Take a copy of the modulus --- */
196
f4535c64 197 if (!MP_POSP(m) || !MP_ODDP(m))
198 return (-1);
d3409d5e 199 mm->m = MP_COPY(m);
200
f5f35081 201 /* --- Determine %$R^2$% --- */
d3409d5e 202
f5f35081 203 mm->n = n;
204 MPX_ZERO(r2->v, r2->vl - 1);
205 r2->vl[-1] = 1;
d3409d5e 206
f5f35081 207 /* --- Find the magic value @mi@ --- */
208
209 mp_build(&r, r2->v + n, r2->vl);
b817bfc6 210 mm->mi = mp_modinv(MP_NEW, m, &r);
f5f35081 211 mm->mi = mp_sub(mm->mi, &r, mm->mi);
362c3d18 212 MP_ENSURE(mm->mi, n);
d3409d5e 213
214 /* --- Discover the values %$R \bmod m$% and %$R^2 \bmod m$% --- */
215
f5f35081 216 mm->r2 = MP_NEW;
217 mp_div(0, &mm->r2, r2, m);
218 mm->r = mpmont_reduce(mm, MP_NEW, mm->r2);
219 MP_DROP(r2);
f4535c64 220 return (0);
d3409d5e 221}
222
52e4b041 223#endif
224
d3409d5e 225/* --- @mpmont_destroy@ --- *
226 *
227 * Arguments: @mpmont *mm@ = pointer to a Montgomery reduction context
228 *
229 * Returns: ---
230 *
231 * Use: Disposes of a context when it's no longer of any use to
232 * anyone.
233 */
234
235void mpmont_destroy(mpmont *mm)
236{
237 MP_DROP(mm->m);
238 MP_DROP(mm->r);
239 MP_DROP(mm->r2);
f5f35081 240 MP_DROP(mm->mi);
d3409d5e 241}
242
243/* --- @mpmont_reduce@ --- *
244 *
245 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
246 * @mp *d@ = destination
ef5f4810 247 * @mp *a@ = source, assumed positive
d3409d5e 248 *
249 * Returns: Result, %$a R^{-1} \bmod m$%.
250 */
251
52e4b041 252#ifdef MPMONT_DISABLE
253
ef5f4810 254mp *mpmont_reduce(mpmont *mm, mp *d, mp *a)
52e4b041 255{
256 mp_div(0, &d, a, mm->m);
257 return (d);
258}
259
260#else
261
ef5f4810 262mp *mpmont_reduce(mpmont *mm, mp *d, mp *a)
d3409d5e 263{
f5f35081 264 size_t n = mm->n;
265
266 /* --- Check for serious Karatsuba reduction --- */
267
52cdaca9 268 if (n > MPK_THRESH * 3) {
f5f35081 269 mp al;
270 mpw *vl;
271 mp *u;
272
0e70bd46
MW
273 if (MP_LEN(a) >= n) vl = a->v + n;
274 else vl = a->vl;
f5f35081 275 mp_build(&al, a->v, vl);
276 u = mp_mul(MP_NEW, &al, mm->mi);
0e70bd46 277 if (MP_LEN(u) > n) u->vl = u->v + n;
f5f35081 278 u = mp_mul(u, u, mm->m);
279 d = mp_add(d, a, u);
366e9701 280 MP_ENSURE(d, n);
f5f35081 281 mp_drop(u);
282 }
d3409d5e 283
f5f35081 284 /* --- Otherwise do it the hard way --- */
d3409d5e 285
d3409d5e 286 else {
d34decd2 287 a = MP_COPY(a);
0e70bd46 288 if (d) MP_DROP(d);
d34decd2 289 d = a;
0e70bd46
MW
290 MP_DEST(d, 2*mm->n + 1, a->f);
291 redccore(d->v, d->vl, mm->m->v, mm->n, mm->mi->v);
d3409d5e 292 }
293
f5f35081 294 /* --- Wrap everything up --- */
d3409d5e 295
0e70bd46 296 finish(mm, d);
d3409d5e 297 return (d);
298}
299
52e4b041 300#endif
301
d3409d5e 302/* --- @mpmont_mul@ --- *
303 *
304 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
305 * @mp *d@ = destination
ef5f4810 306 * @mp *a, *b@ = sources, assumed positive
d3409d5e 307 *
308 * Returns: Result, %$a b R^{-1} \bmod m$%.
309 */
310
52e4b041 311#ifdef MPMONT_DISABLE
312
ef5f4810 313mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b)
52e4b041 314{
315 d = mp_mul(d, a, b);
316 mp_div(0, &d, d, mm->m);
317 return (d);
318}
319
320#else
321
ef5f4810 322mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b)
d3409d5e 323{
52cdaca9 324 if (mm->n > MPK_THRESH * 3) {
ef5f4810 325 d = mp_mul(d, a, b);
326 d = mpmont_reduce(mm, d, d);
327 } else {
ef5f4810 328 a = MP_COPY(a);
329 b = MP_COPY(b);
0e70bd46
MW
330 MP_DEST(d, 2*mm->n + 1, a->f | b->f | MP_UNDEF);
331 mulcore(d->v, d->vl, a->v, a->vl, b->v, b->vl,
332 mm->m->v, mm->n, mm->mi->v);
333 d->f = ((a->f | b->f) & MP_BURN) | ((a->f ^ b->f) & MP_NEG);
334 finish(mm, d);
335 MP_DROP(a); MP_DROP(b);
d3409d5e 336 }
337
d3409d5e 338 return (d);
339}
340
52e4b041 341#endif
342
d3409d5e 343/*----- Test rig ----------------------------------------------------------*/
344
345#ifdef TEST_RIG
346
347static int tcreate(dstr *v)
348{
349 mp *m = *(mp **)v[0].buf;
350 mp *mi = *(mp **)v[1].buf;
351 mp *r = *(mp **)v[2].buf;
352 mp *r2 = *(mp **)v[3].buf;
353
354 mpmont mm;
355 int ok = 1;
356
357 mpmont_create(&mm, m);
358
f5f35081 359 if (mm.mi->v[0] != mi->v[0]) {
d3409d5e 360 fprintf(stderr, "\n*** bad mi: found %lu, expected %lu",
f5f35081 361 (unsigned long)mm.mi->v[0], (unsigned long)mi->v[0]);
d3409d5e 362 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
363 fputc('\n', stderr);
364 ok = 0;
365 }
366
032099d1 367 if (!MP_EQ(mm.r, r)) {
d3409d5e 368 fputs("\n*** bad r", stderr);
369 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
370 fputs("\nexpected ", stderr); mp_writefile(r, stderr, 10);
45c0fd36 371 fputs("\n found ", stderr); mp_writefile(mm.r, stderr, 10);
d3409d5e 372 fputc('\n', stderr);
373 ok = 0;
374 }
375
032099d1 376 if (!MP_EQ(mm.r2, r2)) {
d3409d5e 377 fputs("\n*** bad r2", stderr);
378 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
379 fputs("\nexpected ", stderr); mp_writefile(r2, stderr, 10);
45c0fd36 380 fputs("\n found ", stderr); mp_writefile(mm.r2, stderr, 10);
d3409d5e 381 fputc('\n', stderr);
382 ok = 0;
383 }
384
385 MP_DROP(m);
386 MP_DROP(mi);
387 MP_DROP(r);
388 MP_DROP(r2);
389 mpmont_destroy(&mm);
ef5f4810 390 assert(mparena_count(MPARENA_GLOBAL) == 0);
d3409d5e 391 return (ok);
392}
393
394static int tmul(dstr *v)
395{
396 mp *m = *(mp **)v[0].buf;
397 mp *a = *(mp **)v[1].buf;
398 mp *b = *(mp **)v[2].buf;
399 mp *r = *(mp **)v[3].buf;
d3409d5e 400 int ok = 1;
401
402 mpmont mm;
403 mpmont_create(&mm, m);
404
405 {
79a34029 406 mp *qr = mp_mul(MP_NEW, a, b);
407 mp_div(0, &qr, qr, m);
408
032099d1 409 if (!MP_EQ(qr, r)) {
79a34029 410 fputs("\n*** classical modmul failed", stderr);
411 fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
412 fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
413 fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
414 fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
415 fputs("\nqr = ", stderr); mp_writefile(qr, stderr, 10);
416 fputc('\n', stderr);
417 ok = 0;
418 }
419
420 mp_drop(qr);
421 }
422
423 {
d3409d5e 424 mp *ar = mpmont_mul(&mm, MP_NEW, a, mm.r2);
425 mp *br = mpmont_mul(&mm, MP_NEW, b, mm.r2);
79a34029 426 mp *mr = mpmont_mul(&mm, MP_NEW, ar, br);
d3409d5e 427 mr = mpmont_reduce(&mm, mr, mr);
032099d1 428 if (!MP_EQ(mr, r)) {
79a34029 429 fputs("\n*** montgomery modmul failed", stderr);
430 fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
431 fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
432 fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
433 fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
434 fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10);
435 fputc('\n', stderr);
436 ok = 0;
437 }
d3409d5e 438 MP_DROP(ar); MP_DROP(br);
79a34029 439 mp_drop(mr);
d3409d5e 440 }
441
d3409d5e 442
443 MP_DROP(m);
444 MP_DROP(a);
445 MP_DROP(b);
446 MP_DROP(r);
d3409d5e 447 mpmont_destroy(&mm);
ef5f4810 448 assert(mparena_count(MPARENA_GLOBAL) == 0);
d3409d5e 449 return ok;
450}
451
d3409d5e 452static test_chunk tests[] = {
ef5f4810 453 { "create", tcreate, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
454 { "mul", tmul, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
d3409d5e 455 { 0, 0, { 0 } },
456};
457
458int main(int argc, char *argv[])
459{
460 sub_init();
0f00dc4c 461 test_run(argc, argv, tests, SRCDIR "/t/mpmont");
d3409d5e 462 return (0);
463}
464
465#endif
466
467/*----- That's all, folks -------------------------------------------------*/