|
2
|
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 } |
