1 /* TomsFastMath, a fast ISO C bignum library.
3 * This project is meant to fill in where LibTomMath
4 * falls short. That is speed ;-)
6 * This project is public domain and free for all purposes.
8 * Tom St Denis, tomstdenis@gmail.com
10 #include <tfm_private.h>
12 #ifdef TFM_TIMING_RESISTANT
14 /* timing resistant montgomery ladder based exptmod
16 Based on work by Marc Joye, Sung-Ming Yen, "The Montgomery Powering Ladder", Cryptographic Hardware and Embedded Systems, CHES 2002
18 static int _fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
22 int err, bitcnt, digidx, y;
24 /* now setup montgomery */
25 if ((err = fp_montgomery_setup (P, &mp)) != FP_OKAY) {
32 /* now we need R mod m */
33 fp_montgomery_calc_normalization (&R[0], P);
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 */
42 fp_mulmod (&R[1], &R[0], P, &R[1]);
44 /* for j = t-1 downto 0 do
45 r_!k = R0*R1; r_k = r_k^2
48 /* set initial mode and bit cnt */
54 /* grab next digit as required */
56 /* if digidx == -1 we are out of digits so break */
60 /* read next digit and reset bitcnt */
61 buf = X->dp[digidx--];
62 bitcnt = (int)DIGIT_BIT;
65 /* grab the next msb from the exponent */
66 y = (fp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
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);
74 fp_montgomery_reduce(&R[0], P, mp);
82 * Some restrictions... x must be positive and < b
84 static int _fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
88 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
90 /* find window size */
91 x = fp_count_bits (X);
96 } else if (x <= 140) {
98 } else if (x <= 450) {
105 memset(M, 0, sizeof(M));
107 /* now setup montgomery */
108 if ((err = fp_montgomery_setup (P, &mp)) != FP_OKAY) {
117 * The M table contains powers of the input base, e.g. M[x] = G^x mod P
119 * The first half of the table is not computed though accept for M[0] and M[1]
122 /* now we need R mod m */
123 fp_montgomery_calc_normalization (&res, P);
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 */
132 fp_mulmod (&M[1], &res, P, &M[1]);
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);
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);
147 /* set initial mode and bit cnt */
151 digidx = X->used - 1;
156 /* grab next digit as required */
158 /* if digidx == -1 we are out of digits so break */
162 /* read next digit and reset bitcnt */
163 buf = X->dp[digidx--];
164 bitcnt = (int)DIGIT_BIT;
167 /* grab the next msb from the exponent */
168 y = (fp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
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
176 if (mode == 0 && y == 0) {
180 /* if the bit is zero and mode == 1 then we square */
181 if (mode == 1 && y == 0) {
183 fp_montgomery_reduce(&res, P, mp);
187 /* else we add it to the window */
188 bitbuf |= (y << (winsize - ++bitcpy));
191 if (bitcpy == winsize) {
192 /* ok window is filled so square as required and multiply */
194 for (x = 0; x < winsize; x++) {
196 fp_montgomery_reduce(&res, P, mp);
200 fp_mul(&res, &M[bitbuf], &res);
201 fp_montgomery_reduce(&res, P, mp);
203 /* empty window and reset */
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++) {
215 fp_montgomery_reduce(&res, P, mp);
217 /* get next bit of the window */
219 if ((bitbuf & (1 << winsize)) != 0) {
221 fp_mul(&res, &M[1], &res);
222 fp_montgomery_reduce(&res, P, mp);
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
233 fp_montgomery_reduce(&res, P, mp);
235 /* swap res with Y */
243 int fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
249 /* prevent overflows */
250 if (P->used > (FP_SIZE/2)) {
256 if (X->sign == FP_NEG) {
257 /* yes, copy G and invmod it */
259 if ((err = fp_invmod(&tmp, P, &tmp)) != FP_OKAY) {
263 err = _fp_exptmod(&tmp, X, P, Y);
269 /* Positive exponent so just exptmod */
270 return _fp_exptmod(G, X, P, Y);