Mercurial > dropbear

comparison bn_mp_exptmod_fast.c @ 2:86e0b50a9b58 libtommath-orig ltm-0.30-orig

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ltm 0.30 orig import
author Matt Johnston <matt@ucc.asn.au>
date Mon, 31 May 2004 18:25:22 +0000
parents
children d29b64170cf0
comparison
equal deleted inserted replaced
-1:000000000000 2:86e0b50a9b58
1 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2 *
3 * LibTomMath is a library that provides multiple-precision
4 * integer arithmetic as well as number theoretic functionality.
5 *
6 * The library was designed directly after the MPI library by
7 * Michael Fromberger but has been written from scratch with
8 * additional optimizations in place.
9 *
10 * The library is free for all purposes without any express
11 * guarantee it works.
12 *
13 * Tom St Denis, [email protected], http://math.libtomcrypt.org
14 */
15 #include <tommath.h>
16
17 /* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
18 *
19 * Uses a left-to-right k-ary sliding window to compute the modular exponentiation.
20 * The value of k changes based on the size of the exponent.
21 *
22 * Uses Montgomery or Diminished Radix reduction [whichever appropriate]
23 */
24
25 #ifdef MP_LOW_MEM
26 #define TAB_SIZE 32
27 #else
28 #define TAB_SIZE 256
29 #endif
30
31 int
32 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 /* now setup montgomery */
88 if ((err = mp_montgomery_setup (P, &mp)) != MP_OKAY) {
89 goto __M;
90 }
91
92 /* automatically pick the comba one if available (saves quite a few calls/ifs) */
93 if (((P->used * 2 + 1) < MP_WARRAY) &&
94 P->used < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
95 redux = fast_mp_montgomery_reduce;
96 } else {
97 /* use slower baseline Montgomery method */
98 redux = mp_montgomery_reduce;
99 }
100 } else if (redmode == 1) {
101 /* setup DR reduction for moduli of the form B**k - b */
102 mp_dr_setup(P, &mp);
103 redux = mp_dr_reduce;
104 } else {
105 /* setup DR reduction for moduli of the form 2**k - b */
106 if ((err = mp_reduce_2k_setup(P, &mp)) != MP_OKAY) {
107 goto __M;
108 }
109 redux = mp_reduce_2k;
110 }
111
112 /* setup result */
113 if ((err = mp_init (&res)) != MP_OKAY) {
114 goto __M;
115 }
116
117 /* create M table
118 *
119 * The M table contains powers of the input base, e.g. M[x] = G^x mod P
120 *
121 * The first half of the table is not computed though accept for M[0] and M[1]
122 */
123
124 if (redmode == 0) {
125 /* now we need R mod m */
126 if ((err = mp_montgomery_calc_normalization (&res, P)) != MP_OKAY) {
127 goto __RES;
128 }
129
130 /* now set M[1] to G * R mod m */
131 if ((err = mp_mulmod (G, &res, P, &M[1])) != MP_OKAY) {
132 goto __RES;
133 }
134 } else {
135 mp_set(&res, 1);
136 if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {
137 goto __RES;
138 }
139 }
140
141 /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
142 if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
143 goto __RES;
144 }
145
146 for (x = 0; x < (winsize - 1); x++) {
147 if ((err = mp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) {
148 goto __RES;
149 }
150 if ((err = redux (&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) {
151 goto __RES;
152 }
153 }
154
155 /* create upper table */
156 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
157 if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
158 goto __RES;
159 }
160 if ((err = redux (&M[x], P, mp)) != MP_OKAY) {
161 goto __RES;
162 }
163 }
164
165 /* set initial mode and bit cnt */
166 mode = 0;
167 bitcnt = 1;
168 buf = 0;
169 digidx = X->used - 1;
170 bitcpy = 0;
171 bitbuf = 0;
172
173 for (;;) {
174 /* grab next digit as required */
175 if (--bitcnt == 0) {
176 /* if digidx == -1 we are out of digits so break */
177 if (digidx == -1) {
178 break;
179 }
180 /* read next digit and reset bitcnt */
181 buf = X->dp[digidx--];
182 bitcnt = (int)DIGIT_BIT;
183 }
184
185 /* grab the next msb from the exponent */
186 y = (mp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
187 buf <<= (mp_digit)1;
188
189 /* if the bit is zero and mode == 0 then we ignore it
190 * These represent the leading zero bits before the first 1 bit
191 * in the exponent. Technically this opt is not required but it
192 * does lower the # of trivial squaring/reductions used
193 */
194 if (mode == 0 && y == 0) {
195 continue;
196 }
197
198 /* if the bit is zero and mode == 1 then we square */
199 if (mode == 1 && y == 0) {
200 if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
201 goto __RES;
202 }
203 if ((err = redux (&res, P, mp)) != MP_OKAY) {
204 goto __RES;
205 }
206 continue;
207 }
208
209 /* else we add it to the window */
210 bitbuf |= (y << (winsize - ++bitcpy));
211 mode = 2;
212
213 if (bitcpy == winsize) {
214 /* ok window is filled so square as required and multiply */
215 /* square first */
216 for (x = 0; x < winsize; x++) {
217 if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
218 goto __RES;
219 }
220 if ((err = redux (&res, P, mp)) != MP_OKAY) {
221 goto __RES;
222 }
223 }
224
225 /* then multiply */
226 if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
227 goto __RES;
228 }
229 if ((err = redux (&res, P, mp)) != MP_OKAY) {
230 goto __RES;
231 }
232
233 /* empty window and reset */
234 bitcpy = 0;
235 bitbuf = 0;
236 mode = 1;
237 }
238 }
239
240 /* if bits remain then square/multiply */
241 if (mode == 2 && bitcpy > 0) {
242 /* square then multiply if the bit is set */
243 for (x = 0; x < bitcpy; x++) {
244 if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
245 goto __RES;
246 }
247 if ((err = redux (&res, P, mp)) != MP_OKAY) {
248 goto __RES;
249 }
250
251 /* get next bit of the window */
252 bitbuf <<= 1;
253 if ((bitbuf & (1 << winsize)) != 0) {
254 /* then multiply */
255 if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
256 goto __RES;
257 }
258 if ((err = redux (&res, P, mp)) != MP_OKAY) {
259 goto __RES;
260 }
261 }
262 }
263 }
264
265 if (redmode == 0) {
266 /* fixup result if Montgomery reduction is used
267 * recall that any value in a Montgomery system is
268 * actually multiplied by R mod n. So we have
269 * to reduce one more time to cancel out the factor
270 * of R.
271 */
272 if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) {
273 goto __RES;
274 }
275 }
276
277 /* swap res with Y */
278 mp_exch (&res, Y);
279 err = MP_OKAY;
280 __RES:mp_clear (&res);
281 __M:
282 mp_clear(&M[1]);
283 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
284 mp_clear (&M[x]);
285 }
286 return err;
287 }