3 * $Id: mpmont.c,v 1.3 1999/11/21 11:35:10 mdw Exp $
7 * (c) 1999 Straylight/Edgeware
10 /*----- Licensing notice --------------------------------------------------*
12 * This file is part of Catacomb.
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.
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.
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,
30 /*----- Revision history --------------------------------------------------*
33 * Revision 1.3 1999/11/21 11:35:10 mdw
34 * Performance improvement: use @mp_sqr@ and @mpmont_reduce@ instead of
35 * @mpmont_mul@ for squaring in exponentiation.
37 * Revision 1.2 1999/11/19 13:17:26 mdw
38 * Add extra interface to exponentiation which returns a Montgomerized
41 * Revision 1.1 1999/11/17 18:02:16 mdw
42 * New multiprecision integer arithmetic suite.
46 /*----- Header files ------------------------------------------------------*/
51 /*----- Main code ---------------------------------------------------------*/
53 /* --- @mpmont_create@ --- *
55 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
56 * @mp *m@ = modulus to use
60 * Use: Initializes a Montgomery reduction context ready for use.
63 void mpmont_create(mpmont
*mm
, mp
*m
)
65 /* --- Take a copy of the modulus --- */
70 /* --- Find the magic value @mi@ --- *
72 * This is a slightly grungy way of solving the problem, but it does work.
81 mp_build(&a
, av
, av
+ 2);
82 mp_build(&b
, m
->v
, m
->v
+ 1);
83 mp_gcd(0, 0, &i
, &a
, &b
);
91 /* --- Discover the values %$R \bmod m$% and %$R^2 \bmod m$% --- */
95 mp
*r
= mp_create(l
+ 1);
97 mm
->shift
= l
* MPW_BITS
;
98 MPX_ZERO(r
->v
, r
->vl
- 1);
100 mm
->r
= mm
->r2
= MP_NEW
;
101 mp_div(0, &mm
->r
, r
, m
);
102 r
= mp_sqr(r
, mm
->r
);
103 mp_div(0, &mm
->r2
, r
, m
);
108 /* --- @mpmont_destroy@ --- *
110 * Arguments: @mpmont *mm@ = pointer to a Montgomery reduction context
114 * Use: Disposes of a context when it's no longer of any use to
118 void mpmont_destroy(mpmont
*mm
)
125 /* --- @mpmont_reduce@ --- *
127 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
128 * @mp *d@ = destination
129 * @const mp *a@ = source, assumed positive
131 * Returns: Result, %$a R^{-1} \bmod m$%.
134 mp
*mpmont_reduce(mpmont
*mm
, mp
*d
, const mp
*a
)
140 /* --- Initial conditioning of the arguments --- */
145 MP_MODIFY(d
, 2 * n
+ 1);
147 MP_MODIFY(d
, 2 * n
+ 1);
148 memcpy(d
->v
, a
->v
, MPWS(MP_LEN(a
)));
149 memset(d
->v
+ MP_LEN(a
), 0, MPWS(MP_LEN(d
) - MP_LEN(a
)));
152 dv
= d
->v
; dvl
= d
->vl
;
153 mv
= mm
->m
->v
; mvl
= mm
->m
->vl
;
155 /* --- Let's go to work --- */
158 mpw u
= MPW(*dv
* mm
->mi
);
159 MPX_UMLAN(dv
, dvl
, mv
, mvl
, u
);
165 memmove(d
->v
, dv
, MPWS(dvl
- dv
));
168 d
->f
= a
->f
& MP_BURN
;
169 if (MP_CMP(d
, >=, mm
->m
))
170 d
= mp_sub(d
, d
, mm
->m
);
174 /* --- @mpmont_mul@ --- *
176 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
177 * @mp *d@ = destination
178 * @const mp *a, *b@ = sources, assumed positive
180 * Returns: Result, %$a b R^{-1} \bmod m$%.
183 mp
*mpmont_mul(mpmont
*mm
, mp
*d
, const mp
*a
, const mp
*b
)
192 /* --- Initial conditioning of the arguments --- */
194 if (MP_LEN(a
) > MP_LEN(b
)) {
195 const mp
*t
= a
; a
= b
; b
= t
;
199 MP_MODIFY(d
, 2 * n
+ 1);
200 dv
= d
->v
; dvl
= d
->vl
;
202 av
= a
->v
; avl
= a
->vl
;
203 bv
= b
->v
; bvl
= b
->vl
;
204 mv
= mm
->m
->v
; mvl
= mm
->m
->vl
;
207 /* --- Montgomery multiplication phase --- */
210 while (i
< n
&& av
< avl
) {
212 mpw u
= MPW((*dv
+ x
* y
) * mm
->mi
);
213 MPX_UMLAN(dv
, dvl
, bv
, bvl
, x
);
214 MPX_UMLAN(dv
, dvl
, mv
, mvl
, u
);
219 /* --- Simpler Montgomery reduction phase --- */
222 mpw u
= MPW(*dv
* mm
->mi
);
223 MPX_UMLAN(dv
, dvl
, mv
, mvl
, u
);
230 memmove(d
->v
, dv
, MPWS(dvl
- dv
));
233 d
->f
= (a
->f
| b
->f
) & MP_BURN
;
234 if (MP_CMP(d
, >=, mm
->m
))
235 d
= mp_sub(d
, d
, mm
->m
);
239 /* --- @mpmont_expr@ --- *
241 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
242 * @const mp *a@ = base
243 * @const mp *e@ = exponent
245 * Returns: Result, %$a^e R \bmod m$%.
248 mp
*mpmont_expr(mpmont
*mm
, const mp
*a
, const mp
*e
)
251 mp
*ar
= mpmont_mul(mm
, MP_NEW
, a
, mm
->r2
);
252 mp
*d
= MP_COPY(mm
->r
);
263 dd
= mp_sqr(spare
, ar
);
264 dd
= mpmont_reduce(mm
, dd
, dd
);
268 dd
= mpmont_mul(mm
, spare
, d
, ar
);
282 /* --- @mpmont_exp@ --- *
284 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
285 * @const mp *a@ = base
286 * @const mp *e@ = exponent
288 * Returns: Result, %$a^e \bmod m$%.
291 mp
*mpmont_exp(mpmont
*mm
, const mp
*a
, const mp
*e
)
293 mp
*d
= mpmont_expr(mm
, a
, e
);
294 d
= mpmont_reduce(mm
, d
, d
);
298 /*----- Test rig ----------------------------------------------------------*/
302 static int tcreate(dstr
*v
)
304 mp
*m
= *(mp
**)v
[0].buf
;
305 mp
*mi
= *(mp
**)v
[1].buf
;
306 mp
*r
= *(mp
**)v
[2].buf
;
307 mp
*r2
= *(mp
**)v
[3].buf
;
312 mpmont_create(&mm
, m
);
314 if (mm
.mi
!= mi
->v
[0]) {
315 fprintf(stderr
, "\n*** bad mi: found %lu, expected %lu",
316 (unsigned long)mm
.mi
, (unsigned long)mi
->v
[0]);
317 fputs("\nm = ", stderr
); mp_writefile(m
, stderr
, 10);
322 if (MP_CMP(mm
.r
, !=, r
)) {
323 fputs("\n*** bad r", stderr
);
324 fputs("\nm = ", stderr
); mp_writefile(m
, stderr
, 10);
325 fputs("\nexpected ", stderr
); mp_writefile(r
, stderr
, 10);
326 fputs("\n found ", stderr
); mp_writefile(mm
.r
, stderr
, 10);
331 if (MP_CMP(mm
.r2
, !=, r2
)) {
332 fputs("\n*** bad r2", stderr
);
333 fputs("\nm = ", stderr
); mp_writefile(m
, stderr
, 10);
334 fputs("\nexpected ", stderr
); mp_writefile(r2
, stderr
, 10);
335 fputs("\n found ", stderr
); mp_writefile(mm
.r2
, stderr
, 10);
348 static int tmul(dstr
*v
)
350 mp
*m
= *(mp
**)v
[0].buf
;
351 mp
*a
= *(mp
**)v
[1].buf
;
352 mp
*b
= *(mp
**)v
[2].buf
;
353 mp
*r
= *(mp
**)v
[3].buf
;
357 mpmont_create(&mm
, m
);
360 mp
*qr
= mp_mul(MP_NEW
, a
, b
);
361 mp_div(0, &qr
, qr
, m
);
363 if (MP_CMP(qr
, !=, r
)) {
364 fputs("\n*** classical modmul failed", stderr
);
365 fputs("\n m = ", stderr
); mp_writefile(m
, stderr
, 10);
366 fputs("\n a = ", stderr
); mp_writefile(a
, stderr
, 10);
367 fputs("\n b = ", stderr
); mp_writefile(b
, stderr
, 10);
368 fputs("\n r = ", stderr
); mp_writefile(r
, stderr
, 10);
369 fputs("\nqr = ", stderr
); mp_writefile(qr
, stderr
, 10);
378 mp
*ar
= mpmont_mul(&mm
, MP_NEW
, a
, mm
.r2
);
379 mp
*br
= mpmont_mul(&mm
, MP_NEW
, b
, mm
.r2
);
380 mp
*mr
= mpmont_mul(&mm
, MP_NEW
, ar
, br
);
381 mr
= mpmont_reduce(&mm
, mr
, mr
);
382 if (MP_CMP(mr
, !=, r
)) {
383 fputs("\n*** montgomery modmul failed", stderr
);
384 fputs("\n m = ", stderr
); mp_writefile(m
, stderr
, 10);
385 fputs("\n a = ", stderr
); mp_writefile(a
, stderr
, 10);
386 fputs("\n b = ", stderr
); mp_writefile(b
, stderr
, 10);
387 fputs("\n r = ", stderr
); mp_writefile(r
, stderr
, 10);
388 fputs("\nmr = ", stderr
); mp_writefile(mr
, stderr
, 10);
392 MP_DROP(ar
); MP_DROP(br
);
405 static int texp(dstr
*v
)
407 mp
*m
= *(mp
**)v
[0].buf
;
408 mp
*a
= *(mp
**)v
[1].buf
;
409 mp
*b
= *(mp
**)v
[2].buf
;
410 mp
*r
= *(mp
**)v
[3].buf
;
415 mpmont_create(&mm
, m
);
417 mr
= mpmont_exp(&mm
, a
, b
);
419 if (MP_CMP(mr
, !=, r
)) {
420 fputs("\n*** montgomery modexp failed", stderr
);
421 fputs("\n m = ", stderr
); mp_writefile(m
, stderr
, 10);
422 fputs("\n a = ", stderr
); mp_writefile(a
, stderr
, 10);
423 fputs("\n e = ", stderr
); mp_writefile(b
, stderr
, 10);
424 fputs("\n r = ", stderr
); mp_writefile(r
, stderr
, 10);
425 fputs("\nmr = ", stderr
); mp_writefile(mr
, stderr
, 10);
440 static test_chunk tests
[] = {
441 { "create", tcreate
, { &type_mp
, &type_mp
, &type_mp
, &type_mp
} },
442 { "mul", tmul
, { &type_mp
, &type_mp
, &type_mp
, &type_mp
} },
443 { "exp", texp
, { &type_mp
, &type_mp
, &type_mp
, &type_mp
} },
447 int main(int argc
, char *argv
[])
450 test_run(argc
, argv
, tests
, SRCDIR
"/tests/mpmont");
456 /*----- That's all, folks -------------------------------------------------*/