1 /* TomsFastMath, a fast ISO C bignum library.
2 *
3 * This project is meant to fill in where LibTomMath
4 * falls short. That is speed ;-)
5 *
6 * This project is public domain and free for all purposes.
7 *
8 * Tom St Denis, tomstdenis@gmail.com
9 */
10 #include "bignum_fast.h"
11
12 #ifdef TFM_TIMING_RESISTANT
13
14 /* timing resistant montgomery ladder based exptmod
15
16 Based on work by Marc Joye, Sung-Ming Yen, "The Montgomery Powering Ladder", Cryptographic Hardware and Embedded Systems, CHES 2002
17 */
_fp_exptmod(fp_int * G,fp_int * X,fp_int * P,fp_int * Y)18 static int _fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
19 {
20 fp_int R[2];
21 fp_digit buf, mp;
22 int err, bitcnt, digidx, y;
23
24 /* now setup montgomery */
25 if ((err = fp_montgomery_setup (P, &mp)) != FP_OKAY) {
26 return err;
27 }
28
29 fp_init(&R[0]);
30 fp_init(&R[1]);
31
32 /* now we need R mod m */
33 fp_montgomery_calc_normalization (&R[0], P);
34
35 /* now set R[0][1] to G * R mod m */
36 if (fp_cmp_mag(P, G) != FP_GT) {
37 /* G > P so we reduce it first */
38 fp_mod(G, P, &R[1]);
39 } else {
40 fp_copy(G, &R[1]);
41 }
42 fp_mulmod (&R[1], &R[0], P, &R[1]);
43
44 /* for j = t-1 downto 0 do
45 r_!k = R0*R1; r_k = r_k^2
46 */
47
48 /* set initial mode and bit cnt */
49 bitcnt = 1;
50 buf = 0;
51 digidx = X->used - 1;
52
53 for (;;) {
54 /* grab next digit as required */
55 if (--bitcnt == 0) {
56 /* if digidx == -1 we are out of digits so break */
57 if (digidx == -1) {
58 break;
59 }
60 /* read next digit and reset bitcnt */
61 buf = X->dp[digidx--];
62 bitcnt = (int)DIGIT_BIT;
63 }
64
65 /* grab the next msb from the exponent */
66 y = (fp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
67 buf <<= (fp_digit)1;
68
69 /* do ops */
70 fp_mul(&R[0], &R[1], &R[y^1]); fp_montgomery_reduce(&R[y^1], P, mp);
71 fp_sqr(&R[y], &R[y]); fp_montgomery_reduce(&R[y], P, mp);
72 }
73
74 fp_montgomery_reduce(&R[0], P, mp);
75 fp_copy(&R[0], Y);
76 return FP_OKAY;
77 }
78
79 #else
80
81 /* y = g**x (mod b)
82 * Some restrictions... x must be positive and < b
83 */
_fp_exptmod(fp_int * G,fp_int * X,fp_int * P,fp_int * Y)84 static int _fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
85 {
86 fp_int M[64], res;
87 fp_digit buf, mp;
88 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
89
90 /* find window size */
91 x = fp_count_bits (X);
92 if (x <= 21) {
93 winsize = 1;
94 } else if (x <= 36) {
95 winsize = 3;
96 } else if (x <= 140) {
97 winsize = 4;
98 } else if (x <= 450) {
99 winsize = 5;
100 } else {
101 winsize = 6;
102 }
103
104 /* init M array */
105 memset(M, 0, sizeof(M));
106
107 /* now setup montgomery */
108 if ((err = fp_montgomery_setup (P, &mp)) != FP_OKAY) {
109 return err;
110 }
111
112 /* setup result */
113 fp_init(&res);
114
115 /* create M table
116 *
117 * The M table contains powers of the input base, e.g. M[x] = G^x mod P
118 *
119 * The first half of the table is not computed though accept for M[0] and M[1]
120 */
121
122 /* now we need R mod m */
123 fp_montgomery_calc_normalization (&res, P);
124
125 /* now set M[1] to G * R mod m */
126 if (fp_cmp_mag(P, G) != FP_GT) {
127 /* G > P so we reduce it first */
128 fp_mod(G, P, &M[1]);
129 } else {
130 fp_copy(G, &M[1]);
131 }
132 fp_mulmod (&M[1], &res, P, &M[1]);
133
134 /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
135 fp_copy (&M[1], &M[1 << (winsize - 1)]);
136 for (x = 0; x < (winsize - 1); x++) {
137 fp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)]);
138 fp_montgomery_reduce (&M[1 << (winsize - 1)], P, mp);
139 }
140
141 /* create upper table */
142 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
143 fp_mul(&M[x - 1], &M[1], &M[x]);
144 fp_montgomery_reduce(&M[x], P, mp);
145 }
146
147 /* set initial mode and bit cnt */
148 mode = 0;
149 bitcnt = 1;
150 buf = 0;
151 digidx = X->used - 1;
152 bitcpy = 0;
153 bitbuf = 0;
154
155 for (;;) {
156 /* grab next digit as required */
157 if (--bitcnt == 0) {
158 /* if digidx == -1 we are out of digits so break */
159 if (digidx == -1) {
160 break;
161 }
162 /* read next digit and reset bitcnt */
163 buf = X->dp[digidx--];
164 bitcnt = (int)DIGIT_BIT;
165 }
166
167 /* grab the next msb from the exponent */
168 y = (fp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
169 buf <<= (fp_digit)1;
170
171 /* if the bit is zero and mode == 0 then we ignore it
172 * These represent the leading zero bits before the first 1 bit
173 * in the exponent. Technically this opt is not required but it
174 * does lower the # of trivial squaring/reductions used
175 */
176 if (mode == 0 && y == 0) {
177 continue;
178 }
179
180 /* if the bit is zero and mode == 1 then we square */
181 if (mode == 1 && y == 0) {
182 fp_sqr(&res, &res);
183 fp_montgomery_reduce(&res, P, mp);
184 continue;
185 }
186
187 /* else we add it to the window */
188 bitbuf |= (y << (winsize - ++bitcpy));
189 mode = 2;
190
191 if (bitcpy == winsize) {
192 /* ok window is filled so square as required and multiply */
193 /* square first */
194 for (x = 0; x < winsize; x++) {
195 fp_sqr(&res, &res);
196 fp_montgomery_reduce(&res, P, mp);
197 }
198
199 /* then multiply */
200 fp_mul(&res, &M[bitbuf], &res);
201 fp_montgomery_reduce(&res, P, mp);
202
203 /* empty window and reset */
204 bitcpy = 0;
205 bitbuf = 0;
206 mode = 1;
207 }
208 }
209
210 /* if bits remain then square/multiply */
211 if (mode == 2 && bitcpy > 0) {
212 /* square then multiply if the bit is set */
213 for (x = 0; x < bitcpy; x++) {
214 fp_sqr(&res, &res);
215 fp_montgomery_reduce(&res, P, mp);
216
217 /* get next bit of the window */
218 bitbuf <<= 1;
219 if ((bitbuf & (1 << winsize)) != 0) {
220 /* then multiply */
221 fp_mul(&res, &M[1], &res);
222 fp_montgomery_reduce(&res, P, mp);
223 }
224 }
225 }
226
227 /* fixup result if Montgomery reduction is used
228 * recall that any value in a Montgomery system is
229 * actually multiplied by R mod n. So we have
230 * to reduce one more time to cancel out the factor
231 * of R.
232 */
233 fp_montgomery_reduce(&res, P, mp);
234
235 /* swap res with Y */
236 fp_copy (&res, Y);
237 return FP_OKAY;
238 }
239
240 #endif
241
242
fp_exptmod(fp_int * G,fp_int * X,fp_int * P,fp_int * Y)243 int fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
244 {
245 fp_int tmp;
246 int err;
247
248 #ifdef TFM_CHECK
249 /* prevent overflows */
250 if (P->used > (FP_SIZE/2)) {
251 return FP_VAL;
252 }
253 #endif
254
255 /* is X negative? */
256 if (X->sign == FP_NEG) {
257 /* yes, copy G and invmod it */
258 fp_copy(G, &tmp);
259 if ((err = fp_invmod(&tmp, P, &tmp)) != FP_OKAY) {
260 return err;
261 }
262 X->sign = FP_ZPOS;
263 err = _fp_exptmod(&tmp, X, P, Y);
264 if (X != Y) {
265 X->sign = FP_NEG;
266 }
267 return err;
268 } else {
269 /* Positive exponent so just exptmod */
270 return _fp_exptmod(G, X, P, Y);
271 }
272 }
273
274 /* $Source: /cvs/libtom/tomsfastmath/src/exptmod/fp_exptmod.c,v $ */
275 /* $Revision: 1.1 $ */
276 /* $Date: 2006/12/31 21:25:53 $ */
277