|
282
|
1 #include <tommath.h> |
|
|
2 #ifdef BN_S_MP_EXPTMOD_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 * |
|
386
|
15 * Tom St Denis, [email protected], http://math.libtomcrypt.com |
|
282
|
16 */ |
|
|
17 #ifdef MP_LOW_MEM |
|
|
18 #define TAB_SIZE 32 |
|
|
19 #else |
|
|
20 #define TAB_SIZE 256 |
|
|
21 #endif |
|
|
22 |
|
|
23 int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode) |
|
|
24 { |
|
|
25 mp_int M[TAB_SIZE], res, mu; |
|
|
26 mp_digit buf; |
|
|
27 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize; |
|
|
28 int (*redux)(mp_int*,mp_int*,mp_int*); |
|
|
29 |
|
|
30 /* find window size */ |
|
|
31 x = mp_count_bits (X); |
|
|
32 if (x <= 7) { |
|
|
33 winsize = 2; |
|
|
34 } else if (x <= 36) { |
|
|
35 winsize = 3; |
|
|
36 } else if (x <= 140) { |
|
|
37 winsize = 4; |
|
|
38 } else if (x <= 450) { |
|
|
39 winsize = 5; |
|
|
40 } else if (x <= 1303) { |
|
|
41 winsize = 6; |
|
|
42 } else if (x <= 3529) { |
|
|
43 winsize = 7; |
|
|
44 } else { |
|
|
45 winsize = 8; |
|
|
46 } |
|
|
47 |
|
|
48 #ifdef MP_LOW_MEM |
|
|
49 if (winsize > 5) { |
|
|
50 winsize = 5; |
|
|
51 } |
|
|
52 #endif |
|
|
53 |
|
|
54 /* init M array */ |
|
|
55 /* init first cell */ |
|
|
56 if ((err = mp_init(&M[1])) != MP_OKAY) { |
|
|
57 return err; |
|
|
58 } |
|
|
59 |
|
|
60 /* now init the second half of the array */ |
|
|
61 for (x = 1<<(winsize-1); x < (1 << winsize); x++) { |
|
|
62 if ((err = mp_init(&M[x])) != MP_OKAY) { |
|
|
63 for (y = 1<<(winsize-1); y < x; y++) { |
|
|
64 mp_clear (&M[y]); |
|
|
65 } |
|
|
66 mp_clear(&M[1]); |
|
|
67 return err; |
|
|
68 } |
|
|
69 } |
|
|
70 |
|
|
71 /* create mu, used for Barrett reduction */ |
|
|
72 if ((err = mp_init (&mu)) != MP_OKAY) { |
|
|
73 goto LBL_M; |
|
|
74 } |
|
|
75 |
|
|
76 if (redmode == 0) { |
|
|
77 if ((err = mp_reduce_setup (&mu, P)) != MP_OKAY) { |
|
|
78 goto LBL_MU; |
|
|
79 } |
|
|
80 redux = mp_reduce; |
|
|
81 } else { |
|
|
82 if ((err = mp_reduce_2k_setup_l (P, &mu)) != MP_OKAY) { |
|
|
83 goto LBL_MU; |
|
|
84 } |
|
|
85 redux = mp_reduce_2k_l; |
|
|
86 } |
|
|
87 |
|
|
88 /* create M table |
|
|
89 * |
|
|
90 * The M table contains powers of the base, |
|
|
91 * e.g. M[x] = G**x mod P |
|
|
92 * |
|
|
93 * The first half of the table is not |
|
|
94 * computed though accept for M[0] and M[1] |
|
|
95 */ |
|
|
96 if ((err = mp_mod (G, P, &M[1])) != MP_OKAY) { |
|
|
97 goto LBL_MU; |
|
|
98 } |
|
|
99 |
|
|
100 /* compute the value at M[1<<(winsize-1)] by squaring |
|
|
101 * M[1] (winsize-1) times |
|
|
102 */ |
|
|
103 if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) { |
|
|
104 goto LBL_MU; |
|
|
105 } |
|
|
106 |
|
|
107 for (x = 0; x < (winsize - 1); x++) { |
|
|
108 /* square it */ |
|
|
109 if ((err = mp_sqr (&M[1 << (winsize - 1)], |
|
|
110 &M[1 << (winsize - 1)])) != MP_OKAY) { |
|
|
111 goto LBL_MU; |
|
|
112 } |
|
|
113 |
|
|
114 /* reduce modulo P */ |
|
|
115 if ((err = redux (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) { |
|
|
116 goto LBL_MU; |
|
|
117 } |
|
|
118 } |
|
|
119 |
|
|
120 /* create upper table, that is M[x] = M[x-1] * M[1] (mod P) |
|
|
121 * for x = (2**(winsize - 1) + 1) to (2**winsize - 1) |
|
|
122 */ |
|
|
123 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) { |
|
|
124 if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) { |
|
|
125 goto LBL_MU; |
|
|
126 } |
|
|
127 if ((err = redux (&M[x], P, &mu)) != MP_OKAY) { |
|
|
128 goto LBL_MU; |
|
|
129 } |
|
|
130 } |
|
|
131 |
|
|
132 /* setup result */ |
|
|
133 if ((err = mp_init (&res)) != MP_OKAY) { |
|
|
134 goto LBL_MU; |
|
|
135 } |
|
|
136 mp_set (&res, 1); |
|
|
137 |
|
|
138 /* set initial mode and bit cnt */ |
|
|
139 mode = 0; |
|
|
140 bitcnt = 1; |
|
|
141 buf = 0; |
|
|
142 digidx = X->used - 1; |
|
|
143 bitcpy = 0; |
|
|
144 bitbuf = 0; |
|
|
145 |
|
|
146 for (;;) { |
|
|
147 /* grab next digit as required */ |
|
|
148 if (--bitcnt == 0) { |
|
|
149 /* if digidx == -1 we are out of digits */ |
|
|
150 if (digidx == -1) { |
|
|
151 break; |
|
|
152 } |
|
|
153 /* read next digit and reset the bitcnt */ |
|
|
154 buf = X->dp[digidx--]; |
|
|
155 bitcnt = (int) DIGIT_BIT; |
|
|
156 } |
|
|
157 |
|
|
158 /* grab the next msb from the exponent */ |
|
|
159 y = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1; |
|
|
160 buf <<= (mp_digit)1; |
|
|
161 |
|
|
162 /* if the bit is zero and mode == 0 then we ignore it |
|
|
163 * These represent the leading zero bits before the first 1 bit |
|
|
164 * in the exponent. Technically this opt is not required but it |
|
|
165 * does lower the # of trivial squaring/reductions used |
|
|
166 */ |
|
|
167 if (mode == 0 && y == 0) { |
|
|
168 continue; |
|
|
169 } |
|
|
170 |
|
|
171 /* if the bit is zero and mode == 1 then we square */ |
|
|
172 if (mode == 1 && y == 0) { |
|
|
173 if ((err = mp_sqr (&res, &res)) != MP_OKAY) { |
|
|
174 goto LBL_RES; |
|
|
175 } |
|
|
176 if ((err = redux (&res, P, &mu)) != MP_OKAY) { |
|
|
177 goto LBL_RES; |
|
|
178 } |
|
|
179 continue; |
|
|
180 } |
|
|
181 |
|
|
182 /* else we add it to the window */ |
|
|
183 bitbuf |= (y << (winsize - ++bitcpy)); |
|
|
184 mode = 2; |
|
|
185 |
|
|
186 if (bitcpy == winsize) { |
|
|
187 /* ok window is filled so square as required and multiply */ |
|
|
188 /* square first */ |
|
|
189 for (x = 0; x < winsize; x++) { |
|
|
190 if ((err = mp_sqr (&res, &res)) != MP_OKAY) { |
|
|
191 goto LBL_RES; |
|
|
192 } |
|
|
193 if ((err = redux (&res, P, &mu)) != MP_OKAY) { |
|
|
194 goto LBL_RES; |
|
|
195 } |
|
|
196 } |
|
|
197 |
|
|
198 /* then multiply */ |
|
|
199 if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) { |
|
|
200 goto LBL_RES; |
|
|
201 } |
|
|
202 if ((err = redux (&res, P, &mu)) != MP_OKAY) { |
|
|
203 goto LBL_RES; |
|
|
204 } |
|
|
205 |
|
|
206 /* empty window and reset */ |
|
|
207 bitcpy = 0; |
|
|
208 bitbuf = 0; |
|
|
209 mode = 1; |
|
|
210 } |
|
|
211 } |
|
|
212 |
|
|
213 /* if bits remain then square/multiply */ |
|
|
214 if (mode == 2 && bitcpy > 0) { |
|
|
215 /* square then multiply if the bit is set */ |
|
|
216 for (x = 0; x < bitcpy; x++) { |
|
|
217 if ((err = mp_sqr (&res, &res)) != MP_OKAY) { |
|
|
218 goto LBL_RES; |
|
|
219 } |
|
|
220 if ((err = redux (&res, P, &mu)) != MP_OKAY) { |
|
|
221 goto LBL_RES; |
|
|
222 } |
|
|
223 |
|
|
224 bitbuf <<= 1; |
|
|
225 if ((bitbuf & (1 << winsize)) != 0) { |
|
|
226 /* then multiply */ |
|
|
227 if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) { |
|
|
228 goto LBL_RES; |
|
|
229 } |
|
|
230 if ((err = redux (&res, P, &mu)) != MP_OKAY) { |
|
|
231 goto LBL_RES; |
|
|
232 } |
|
|
233 } |
|
|
234 } |
|
|
235 } |
|
|
236 |
|
|
237 mp_exch (&res, Y); |
|
|
238 err = MP_OKAY; |
|
|
239 LBL_RES:mp_clear (&res); |
|
|
240 LBL_MU:mp_clear (&mu); |
|
|
241 LBL_M: |
|
|
242 mp_clear(&M[1]); |
|
|
243 for (x = 1<<(winsize-1); x < (1 << winsize); x++) { |
|
|
244 mp_clear (&M[x]); |
|
|
245 } |
|
|
246 return err; |
|
|
247 } |
|
|
248 #endif |
|
386
|
249 |
|
|
250 /* $Source: /cvs/libtom/libtommath/bn_s_mp_exptmod.c,v $ */ |
|
|
251 /* $Revision: 1.4 $ */ |
|
|
252 /* $Date: 2006/03/31 14:18:44 $ */ |
