Mercurial > dropbear

comparison libtommath/bn_mp_exptmod_fast.c @ 389:5ff8218bcee9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
propagate from branch 'au.asn.ucc.matt.ltm.dropbear' (head 2af95f00ebd5bb7a28b3817db1218442c935388e) to branch 'au.asn.ucc.matt.dropbear' (head ecd779509ef23a8cdf64888904fc9b31d78aa933)
author Matt Johnston <matt@ucc.asn.au>
date Thu, 11 Jan 2007 03:14:55 +0000
parents eed26cff980b
children 4fbf9a7556ed
comparison
equal deleted inserted replaced
388:fb54020f78e1 389:5ff8218bcee9
1 #include <tommath.h>
2 #ifdef BN_MP_EXPTMOD_FAST_C
3 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4 *
5 * LibTomMath is a library that provides multiple-precision
6 * integer arithmetic as well as number theoretic functionality.
7 *
8 * The library was designed directly after the MPI library by
9 * Michael Fromberger but has been written from scratch with
10 * additional optimizations in place.
11 *
12 * The library is free for all purposes without any express
13 * guarantee it works.
14 *
15 * Tom St Denis, [email protected], http://math.libtomcrypt.com
16 */
17
18 /* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
19 *
20 * Uses a left-to-right k-ary sliding window to compute the modular exponentiation.
21 * The value of k changes based on the size of the exponent.
22 *
23 * Uses Montgomery or Diminished Radix reduction [whichever appropriate]
24 */
25
26 #ifdef MP_LOW_MEM
27 #define TAB_SIZE 32
28 #else
29 #define TAB_SIZE 256
30 #endif
31
32 int mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
33 {
34 mp_int M[TAB_SIZE], res;
35 mp_digit buf, mp;
36 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
37
38 /* use a pointer to the reduction algorithm. This allows us to use
39 * one of many reduction algorithms without modding the guts of
40 * the code with if statements everywhere.
41 */
42 int (*redux)(mp_int*,mp_int*,mp_digit);
43
44 /* find window size */
45 x = mp_count_bits (X);
46 if (x <= 7) {
47 winsize = 2;
48 } else if (x <= 36) {
49 winsize = 3;
50 } else if (x <= 140) {
51 winsize = 4;
52 } else if (x <= 450) {
53 winsize = 5;
54 } else if (x <= 1303) {
55 winsize = 6;
56 } else if (x <= 3529) {
57 winsize = 7;
58 } else {
59 winsize = 8;
60 }
61
62 #ifdef MP_LOW_MEM
63 if (winsize > 5) {
64 winsize = 5;
65 }
66 #endif
67
68 /* init M array */
69 /* init first cell */
70 if ((err = mp_init(&M[1])) != MP_OKAY) {
71 return err;
72 }
73
74 /* now init the second half of the array */
75 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
76 if ((err = mp_init(&M[x])) != MP_OKAY) {
77 for (y = 1<<(winsize-1); y < x; y++) {
78 mp_clear (&M[y]);
79 }
80 mp_clear(&M[1]);
81 return err;
82 }
83 }
84
85 /* determine and setup reduction code */
86 if (redmode == 0) {
87 #ifdef BN_MP_MONTGOMERY_SETUP_C
88 /* now setup montgomery */
89 if ((err = mp_montgomery_setup (P, &mp)) != MP_OKAY) {
90 goto LBL_M;
91 }
92 #else
93 err = MP_VAL;
94 goto LBL_M;
95 #endif
96
97 /* automatically pick the comba one if available (saves quite a few calls/ifs) */
98 #ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C
99 if (((P->used * 2 + 1) < MP_WARRAY) &&
100 P->used < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
101 redux = fast_mp_montgomery_reduce;
102 } else
103 #endif
104 {
105 #ifdef BN_MP_MONTGOMERY_REDUCE_C
106 /* use slower baseline Montgomery method */
107 redux = mp_montgomery_reduce;
108 #else
109 err = MP_VAL;
110 goto LBL_M;
111 #endif
112 }
113 } else if (redmode == 1) {
114 #if defined(BN_MP_DR_SETUP_C) && defined(BN_MP_DR_REDUCE_C)
115 /* setup DR reduction for moduli of the form B**k - b */
116 mp_dr_setup(P, &mp);
117 redux = mp_dr_reduce;
118 #else
119 err = MP_VAL;
120 goto LBL_M;
121 #endif
122 } else {
123 #if defined(BN_MP_REDUCE_2K_SETUP_C) && defined(BN_MP_REDUCE_2K_C)
124 /* setup DR reduction for moduli of the form 2**k - b */
125 if ((err = mp_reduce_2k_setup(P, &mp)) != MP_OKAY) {
126 goto LBL_M;
127 }
128 redux = mp_reduce_2k;
129 #else
130 err = MP_VAL;
131 goto LBL_M;
132 #endif
133 }
134
135 /* setup result */
136 if ((err = mp_init (&res)) != MP_OKAY) {
137 goto LBL_M;
138 }
139
140 /* create M table
141 *
142
143 *
144 * The first half of the table is not computed though accept for M[0] and M[1]
145 */
146
147 if (redmode == 0) {
148 #ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C
149 /* now we need R mod m */
150 if ((err = mp_montgomery_calc_normalization (&res, P)) != MP_OKAY) {
151 goto LBL_RES;
152 }
153 #else
154 err = MP_VAL;
155 goto LBL_RES;
156 #endif
157
158 /* now set M[1] to G * R mod m */
159 if ((err = mp_mulmod (G, &res, P, &M[1])) != MP_OKAY) {
160 goto LBL_RES;
161 }
162 } else {
163 mp_set(&res, 1);
164 if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {
165 goto LBL_RES;
166 }
167 }
168
169 /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
170 if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
171 goto LBL_RES;
172 }
173
174 for (x = 0; x < (winsize - 1); x++) {
175 if ((err = mp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) {
176 goto LBL_RES;
177 }
178 if ((err = redux (&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) {
179 goto LBL_RES;
180 }
181 }
182
183 /* create upper table */
184 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
185 if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
186 goto LBL_RES;
187 }
188 if ((err = redux (&M[x], P, mp)) != MP_OKAY) {
189 goto LBL_RES;
190 }
191 }
192
193 /* set initial mode and bit cnt */
194 mode = 0;
195 bitcnt = 1;
196 buf = 0;
197 digidx = X->used - 1;
198 bitcpy = 0;
199 bitbuf = 0;
200
201 for (;;) {
202 /* grab next digit as required */
203 if (--bitcnt == 0) {
204 /* if digidx == -1 we are out of digits so break */
205 if (digidx == -1) {
206 break;
207 }
208 /* read next digit and reset bitcnt */
209 buf = X->dp[digidx--];
210 bitcnt = (int)DIGIT_BIT;
211 }
212
213 /* grab the next msb from the exponent */
214 y = (mp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
215 buf <<= (mp_digit)1;
216
217 /* if the bit is zero and mode == 0 then we ignore it
218 * These represent the leading zero bits before the first 1 bit
219 * in the exponent. Technically this opt is not required but it
220 * does lower the # of trivial squaring/reductions used
221 */
222 if (mode == 0 && y == 0) {
223 continue;
224 }
225
226 /* if the bit is zero and mode == 1 then we square */
227 if (mode == 1 && y == 0) {
228 if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
229 goto LBL_RES;
230 }
231 if ((err = redux (&res, P, mp)) != MP_OKAY) {
232 goto LBL_RES;
233 }
234 continue;
235 }
236
237 /* else we add it to the window */
238 bitbuf |= (y << (winsize - ++bitcpy));
239 mode = 2;
240
241 if (bitcpy == winsize) {
242 /* ok window is filled so square as required and multiply */
243 /* square first */
244 for (x = 0; x < winsize; x++) {
245 if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
246 goto LBL_RES;
247 }
248 if ((err = redux (&res, P, mp)) != MP_OKAY) {
249 goto LBL_RES;
250 }
251 }
252
253 /* then multiply */
254 if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
255 goto LBL_RES;
256 }
257 if ((err = redux (&res, P, mp)) != MP_OKAY) {
258 goto LBL_RES;
259 }
260
261 /* empty window and reset */
262 bitcpy = 0;
263 bitbuf = 0;
264 mode = 1;
265 }
266 }
267
268 /* if bits remain then square/multiply */
269 if (mode == 2 && bitcpy > 0) {
270 /* square then multiply if the bit is set */
271 for (x = 0; x < bitcpy; x++) {
272 if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
273 goto LBL_RES;
274 }
275 if ((err = redux (&res, P, mp)) != MP_OKAY) {
276 goto LBL_RES;
277 }
278
279 /* get next bit of the window */
280 bitbuf <<= 1;
281 if ((bitbuf & (1 << winsize)) != 0) {
282 /* then multiply */
283 if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
284 goto LBL_RES;
285 }
286 if ((err = redux (&res, P, mp)) != MP_OKAY) {
287 goto LBL_RES;
288 }
289 }
290 }
291 }
292
293 if (redmode == 0) {
294 /* fixup result if Montgomery reduction is used
295 * recall that any value in a Montgomery system is
296 * actually multiplied by R mod n. So we have
297 * to reduce one more time to cancel out the factor
298 * of R.
299 */
300 if ((err = redux(&res, P, mp)) != MP_OKAY) {
301 goto LBL_RES;
302 }
303 }
304
305 /* swap res with Y */
306 mp_exch (&res, Y);
307 err = MP_OKAY;
308 LBL_RES:mp_clear (&res);
309 LBL_M:
310 mp_clear(&M[1]);
311 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
312 mp_clear (&M[x]);
313 }
314 return err;
315 }
316 #endif
317
318
319 /* $Source: /cvs/libtom/libtommath/bn_mp_exptmod_fast.c,v $ */
320 /* $Revision: 1.3 $ */
321 /* $Date: 2006/03/31 14:18:44 $ */