3 * $Id: mpbarrett.c,v 1.7 2001/04/19 18:25:26 mdw Exp $
5 * Barrett modular reduction
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 --------------------------------------------------*
32 * $Log: mpbarrett.c,v $
33 * Revision 1.7 2001/04/19 18:25:26 mdw
34 * Use sliding-window exponentiation.
36 * Revision 1.6 2000/10/08 12:03:44 mdw
37 * (mpbarrett_reduce): Cope with negative numbers.
39 * Revision 1.5 2000/07/29 17:04:33 mdw
40 * Change to use left-to-right bitwise exponentiation. This will improve
41 * performance when the base is small.
43 * Revision 1.4 2000/06/17 11:45:09 mdw
44 * Major memory management overhaul. Added arena support. Use the secure
45 * arena for secret integers. Replace and improve the MP management macros
46 * (e.g., replace MP_MODIFY by MP_DEST).
48 * Revision 1.3 1999/12/12 15:08:52 mdw
49 * Don't bother shifting %$q$% in @mpbarrett_reduce@, just skip the least
52 * Revision 1.2 1999/12/11 01:50:56 mdw
53 * Improve initialization slightly.
55 * Revision 1.1 1999/12/10 23:21:59 mdw
56 * Barrett reduction support: works with even moduli.
60 /*----- Header files ------------------------------------------------------*/
63 #include "mpbarrett.h"
65 /*----- Main code ---------------------------------------------------------*/
67 /* --- @mpbarrett_create@ --- *
69 * Arguments: @mpbarrett *mb@ = pointer to Barrett reduction context
70 * @mp *m@ = modulus to work to
75 * Use: Initializes a Barrett reduction context ready for use.
78 void mpbarrett_create(mpbarrett
*mb
, mp
*m
)
82 /* --- Validate the arguments --- */
84 assert(((void)"Barrett modulus must be positive", (m
->f
& MP_NEG
) == 0));
86 /* --- Compute %$\mu$% --- */
91 b
= mp_new(2 * mb
->k
+ 1, 0);
92 MPX_ZERO(b
->v
, b
->vl
- 1);
98 /* --- @mpbarrett_destroy@ --- *
100 * Arguments: @mpbarrett *mb@ = pointer to Barrett reduction context
104 * Use: Destroys a Barrett reduction context releasing any resources
108 void mpbarrett_destroy(mpbarrett
*mb
)
114 /* --- @mpbarrett_reduce@ --- *
116 * Arguments: @mpbarrett *mb@ = pointer to Barrett reduction context
117 * @mp *d@ = destination for result
118 * @mp *m@ = number to reduce
120 * Returns: The residue of @m@ modulo the number in the reduction
123 * Use: Performs an efficient modular reduction.
126 mp
*mpbarrett_reduce(mpbarrett
*mb
, mp
*d
, mp
*m
)
131 /* --- Special case if @m@ is too small --- */
140 /* --- First stage --- */
144 mp_build(&qq
, m
->v
+ (k
- 1), m
->vl
);
145 q
= mp_mul(MP_NEW
, &qq
, mb
->mu
);
146 if (MP_LEN(q
) <= k
) {
154 /* --- Second stage --- */
161 if (MP_LEN(m
) <= k
+ 1)
165 r
= mp_new(k
+ 1, (q
->f
| mb
->m
->f
) & MP_BURN
);
166 mpx_umul(r
->v
, r
->vl
, q
->v
+ k
+ 1, q
->vl
, mb
->m
->v
, mb
->m
->vl
);
167 MP_DEST(d
, k
+ 1, r
->f
);
168 mpx_usub(d
->v
, d
->vl
, m
->v
, mvl
, r
->v
, r
->vl
);
169 d
->f
= (m
->f
| r
->f
) & (MP_BURN
| MP_NEG
);
175 /* --- Final stage --- */
178 while (MPX_UCMP(d
->v
, d
->vl
, >=, mb
->m
->v
, mb
->m
->vl
))
179 mpx_usub(d
->v
, d
->vl
, d
->v
, d
->vl
, mb
->m
->v
, mb
->m
->vl
);
181 /* --- Fix up the sign --- */
184 mpx_usub(d
->v
, d
->vl
, mb
->m
->v
, mb
->m
->vl
, d
->v
, d
->vl
);
192 /* --- @mpbarrett_exp@ --- *
194 * Arguments: @mpbarrett *mb@ = pointer to Barrett reduction context
195 * @mp *d@ = fake destination
199 * Returns: Result, %$a^e \bmod m$%.
203 #define TABSZ (1 << (WINSZ - 1))
205 #define THRESH (((MPW_BITS / WINSZ) << 2) + 1)
207 static mp
*exp_simple(mpbarrett
*mb
, mp
*d
, mp
*a
, mp
*e
)
211 mp
*spare
= (e
->f
& MP_BURN
) ? MP_NEWSEC
: MP_NEW
;
218 while (!MP_RBIT(&sc
))
221 /* --- Do the main body of the work --- */
227 y
= mp_sqr(spare
, x
);
228 y
= mpbarrett_reduce(mb
, y
, y
);
233 mp
*y
= mp_mul(spare
, x
, a
);
234 y
= mpbarrett_reduce(mb
, y
, y
);
246 /* --- Do a final round of squaring --- */
251 y
= mp_sqr(spare
, x
);
252 y
= mpbarrett_reduce(mb
, y
, y
);
266 mp
*mpbarrett_exp(mpbarrett
*mb
, mp
*d
, mp
*a
, mp
*e
)
270 mp
*spare
= (e
->f
& MP_BURN
) ? MP_NEWSEC
: MP_NEW
;
275 /* --- Do we bother? --- */
280 if (MP_LEN(e
) < THRESH
) {
282 return (exp_simple(mb
, d
, a
, e
));
285 /* --- Do the precomputation --- */
287 a2
= mp_sqr(MP_NEW
, a
);
288 a2
= mpbarrett_reduce(mb
, a2
, a2
);
289 tab
= xmalloc(TABSZ
* sizeof(mp
*));
291 for (i
= 1; i
< TABSZ
; i
++) {
292 mp
*x
= mp_mul(MP_NEW
, tab
[i
- 1], a2
);
293 tab
[i
] = mpbarrett_reduce(mb
, x
, x
);
298 /* --- Skip top-end zero bits --- *
300 * If the initial step worked, there must be a set bit somewhere, so keep
301 * stepping until I find it.
305 while (!MP_RBIT(&sc
))
308 /* --- Now for the main work --- */
314 /* --- The next bit is set, so read a window index --- *
316 * Reset @i@ to zero and increment @sq@. Then, until either I read
317 * @WINSZ@ bits or I run out of bits, scan in a bit: if it's clear, bump
318 * the @z@ counter; if it's set, push a set bit into @i@, shift it over
319 * by @z@ bits, bump @sq@ by @z + 1@ and clear @z@. By the end of this
320 * palaver, @i@ is an index to the precomputed value in @tab@.
327 if (l
>= WINSZ
|| !MP_RSTEP(&sc
))
332 i
= ((i
<< 1) | 1) << z
;
338 /* --- Do the squaring --- *
340 * Remember that @sq@ carries over from the zero-skipping stuff below.
345 y
= mp_sqr(spare
, x
);
346 y
= mpbarrett_reduce(mb
, y
, y
);
351 /* --- Do the multiply --- */
353 { mp
*y
= mp_mul(spare
, x
, tab
[i
]); spare
= x
;
354 x
= mpbarrett_reduce(mb
, y
, y
); }
356 /* --- Now grind along through the rest of the bits --- */
368 /* --- Do a final round of squaring --- */
373 y
= mp_sqr(spare
, x
);
374 y
= mpbarrett_reduce(mb
, y
, y
);
381 for (i
= 0; i
< TABSZ
; i
++)
390 /*----- Test rig ----------------------------------------------------------*/
394 static int vmod(dstr
*v
)
396 mp
*x
= *(mp
**)v
[0].buf
;
397 mp
*n
= *(mp
**)v
[1].buf
;
398 mp
*r
= *(mp
**)v
[2].buf
;
403 mpbarrett_create(&mb
, n
);
404 s
= mpbarrett_reduce(&mb
, MP_NEW
, x
);
407 fputs("\n*** barrett reduction failure\n", stderr
);
408 fputs("x = ", stderr
); mp_writefile(x
, stderr
, 10); fputc('\n', stderr
);
409 fputs("n = ", stderr
); mp_writefile(n
, stderr
, 10); fputc('\n', stderr
);
410 fputs("r = ", stderr
); mp_writefile(r
, stderr
, 10); fputc('\n', stderr
);
411 fputs("s = ", stderr
); mp_writefile(s
, stderr
, 10); fputc('\n', stderr
);
415 mpbarrett_destroy(&mb
);
420 assert(mparena_count(MPARENA_GLOBAL
) == 0);
424 static int vexp(dstr
*v
)
426 mp
*m
= *(mp
**)v
[0].buf
;
427 mp
*a
= *(mp
**)v
[1].buf
;
428 mp
*b
= *(mp
**)v
[2].buf
;
429 mp
*r
= *(mp
**)v
[3].buf
;
434 mpbarrett_create(&mb
, m
);
436 mr
= mpbarrett_exp(&mb
, MP_NEW
, a
, b
);
439 fputs("\n*** barrett modexp failed", stderr
);
440 fputs("\n m = ", stderr
); mp_writefile(m
, stderr
, 10);
441 fputs("\n a = ", stderr
); mp_writefile(a
, stderr
, 10);
442 fputs("\n e = ", stderr
); mp_writefile(b
, stderr
, 10);
443 fputs("\n r = ", stderr
); mp_writefile(r
, stderr
, 10);
444 fputs("\nmr = ", stderr
); mp_writefile(mr
, stderr
, 10);
454 mpbarrett_destroy(&mb
);
455 assert(mparena_count(MPARENA_GLOBAL
) == 0);
459 static test_chunk tests
[] = {
460 { "mpbarrett-reduce", vmod
, { &type_mp
, &type_mp
, &type_mp
, 0 } },
461 { "mpbarrett-exp", vexp
, { &type_mp
, &type_mp
, &type_mp
, &type_mp
, 0 } },
465 int main(int argc
, char *argv
[])
468 test_run(argc
, argv
, tests
, SRCDIR
"/tests/mpbarrett");
474 /*----- That's all, folks -------------------------------------------------*/