Line data Source code
1 : #include "tommath_private.h"
2 : #ifdef BN_S_MP_EXPTMOD_FAST_C
3 : /* LibTomMath, multiple-precision integer library -- Tom St Denis */
4 : /* SPDX-License-Identifier: Unlicense */
5 :
6 : /* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
7 : *
8 : * Uses a left-to-right k-ary sliding window to compute the modular exponentiation.
9 : * The value of k changes based on the size of the exponent.
10 : *
11 : * Uses Montgomery or Diminished Radix reduction [whichever appropriate]
12 : */
13 :
14 : #ifdef MP_LOW_MEM
15 : # define TAB_SIZE 32
16 : # define MAX_WINSIZE 5
17 : #else
18 : # define TAB_SIZE 256
19 : # define MAX_WINSIZE 0
20 : #endif
21 :
22 922 : mp_err s_mp_exptmod_fast(const mp_int *G, const mp_int *X, const mp_int *P, mp_int *Y, int redmode)
23 : {
24 32 : mp_int M[TAB_SIZE], res;
25 32 : mp_digit buf, mp;
26 32 : int bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
27 32 : mp_err err;
28 :
29 : /* use a pointer to the reduction algorithm. This allows us to use
30 : * one of many reduction algorithms without modding the guts of
31 : * the code with if statements everywhere.
32 : */
33 32 : mp_err(*redux)(mp_int *x, const mp_int *n, mp_digit rho);
34 :
35 : /* find window size */
36 922 : x = mp_count_bits(X);
37 922 : if (x <= 7) {
38 0 : winsize = 2;
39 922 : } else if (x <= 36) {
40 443 : winsize = 3;
41 463 : } else if (x <= 140) {
42 0 : winsize = 4;
43 463 : } else if (x <= 450) {
44 0 : winsize = 5;
45 463 : } else if (x <= 1303) {
46 243 : winsize = 6;
47 220 : } else if (x <= 3529) {
48 204 : winsize = 7;
49 : } else {
50 0 : winsize = 8;
51 : }
52 :
53 922 : winsize = MAX_WINSIZE ? MP_MIN(MAX_WINSIZE, winsize) : winsize;
54 :
55 : /* init M array */
56 : /* init first cell */
57 922 : if ((err = mp_init_size(&M[1], P->alloc)) != MP_OKAY) {
58 0 : return err;
59 : }
60 :
61 : /* now init the second half of the array */
62 24614 : for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
63 23692 : if ((err = mp_init_size(&M[x], P->alloc)) != MP_OKAY) {
64 0 : for (y = 1<<(winsize-1); y < x; y++) {
65 0 : mp_clear(&M[y]);
66 : }
67 0 : mp_clear(&M[1]);
68 0 : return err;
69 : }
70 : }
71 :
72 : /* determine and setup reduction code */
73 922 : if (redmode == 0) {
74 32 : if (MP_HAS(MP_MONTGOMERY_SETUP)) {
75 : /* now setup montgomery */
76 922 : if ((err = mp_montgomery_setup(P, &mp)) != MP_OKAY) goto LBL_M;
77 : } else {
78 : err = MP_VAL;
79 : goto LBL_M;
80 : }
81 :
82 : /* automatically pick the comba one if available (saves quite a few calls/ifs) */
83 922 : if (MP_HAS(S_MP_MONTGOMERY_REDUCE_FAST) &&
84 922 : (((P->used * 2) + 1) < MP_WARRAY) &&
85 890 : (P->used < MP_MAXFAST)) {
86 890 : redux = s_mp_montgomery_reduce_fast;
87 0 : } else if (MP_HAS(MP_MONTGOMERY_REDUCE)) {
88 : /* use slower baseline Montgomery method */
89 0 : redux = mp_montgomery_reduce;
90 : } else {
91 : err = MP_VAL;
92 : goto LBL_M;
93 : }
94 0 : } else if (redmode == 1) {
95 0 : if (MP_HAS(MP_DR_SETUP) && MP_HAS(MP_DR_REDUCE)) {
96 : /* setup DR reduction for moduli of the form B**k - b */
97 0 : mp_dr_setup(P, &mp);
98 0 : redux = mp_dr_reduce;
99 : } else {
100 : err = MP_VAL;
101 : goto LBL_M;
102 : }
103 0 : } else if (MP_HAS(MP_REDUCE_2K_SETUP) && MP_HAS(MP_REDUCE_2K)) {
104 : /* setup DR reduction for moduli of the form 2**k - b */
105 0 : if ((err = mp_reduce_2k_setup(P, &mp)) != MP_OKAY) goto LBL_M;
106 0 : redux = mp_reduce_2k;
107 : } else {
108 : err = MP_VAL;
109 : goto LBL_M;
110 : }
111 :
112 : /* setup result */
113 922 : if ((err = mp_init_size(&res, P->alloc)) != MP_OKAY) goto LBL_M;
114 :
115 : /* create M table
116 : *
117 :
118 : *
119 : * The first half of the table is not computed though accept for M[0] and M[1]
120 : */
121 :
122 922 : if (redmode == 0) {
123 32 : if (MP_HAS(MP_MONTGOMERY_CALC_NORMALIZATION)) {
124 : /* now we need R mod m */
125 922 : if ((err = mp_montgomery_calc_normalization(&res, P)) != MP_OKAY) goto LBL_RES;
126 :
127 : /* now set M[1] to G * R mod m */
128 922 : if ((err = mp_mulmod(G, &res, P, &M[1])) != MP_OKAY) goto LBL_RES;
129 : } else {
130 : err = MP_VAL;
131 : goto LBL_RES;
132 : }
133 : } else {
134 0 : mp_set(&res, 1uL);
135 0 : if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) goto LBL_RES;
136 : }
137 :
138 : /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
139 922 : if ((err = mp_copy(&M[1], &M[(size_t)1 << (winsize - 1)])) != MP_OKAY) goto LBL_RES;
140 :
141 4375 : for (x = 0; x < (winsize - 1); x++) {
142 3453 : if ((err = mp_sqr(&M[(size_t)1 << (winsize - 1)], &M[(size_t)1 << (winsize - 1)])) != MP_OKAY) goto LBL_RES;
143 3453 : if ((err = redux(&M[(size_t)1 << (winsize - 1)], P, mp)) != MP_OKAY) goto LBL_RES;
144 : }
145 :
146 : /* create upper table */
147 23692 : for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
148 22770 : if ((err = mp_mul(&M[x - 1], &M[1], &M[x])) != MP_OKAY) goto LBL_RES;
149 22770 : if ((err = redux(&M[x], P, mp)) != MP_OKAY) goto LBL_RES;
150 : }
151 :
152 : /* set initial mode and bit cnt */
153 922 : mode = 0;
154 922 : bitcnt = 1;
155 922 : buf = 0;
156 922 : digidx = X->used - 1;
157 922 : bitcpy = 0;
158 922 : bitbuf = 0;
159 :
160 34592 : for (;;) {
161 : /* grab next digit as required */
162 752902 : if (--bitcnt == 0) {
163 : /* if digidx == -1 we are out of digits so break */
164 13455 : if (digidx == -1) {
165 890 : break;
166 : }
167 : /* read next digit and reset bitcnt */
168 12533 : buf = X->dp[digidx--];
169 12533 : bitcnt = (int)MP_DIGIT_BIT;
170 : }
171 :
172 : /* grab the next msb from the exponent */
173 751980 : y = (mp_digit)(buf >> (MP_DIGIT_BIT - 1)) & 1uL;
174 751980 : buf <<= (mp_digit)1;
175 :
176 : /* if the bit is zero and mode == 0 then we ignore it
177 : * These represent the leading zero bits before the first 1 bit
178 : * in the exponent. Technically this opt is not required but it
179 : * does lower the # of trivial squaring/reductions used
180 : */
181 751980 : if ((mode == 0) && (y == 0)) {
182 45308 : continue;
183 : }
184 :
185 : /* if the bit is zero and mode == 1 then we square */
186 706672 : if ((mode == 1) && (y == 0)) {
187 97615 : if ((err = mp_sqr(&res, &res)) != MP_OKAY) goto LBL_RES;
188 97615 : if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
189 97615 : continue;
190 : }
191 :
192 : /* else we add it to the window */
193 609057 : bitbuf |= (y << (winsize - ++bitcpy));
194 609057 : mode = 2;
195 :
196 609057 : if (bitcpy == winsize) {
197 : /* ok window is filled so square as required and multiply */
198 : /* square first */
199 699497 : for (x = 0; x < winsize; x++) {
200 607402 : if ((err = mp_sqr(&res, &res)) != MP_OKAY) goto LBL_RES;
201 607402 : if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
202 : }
203 :
204 : /* then multiply */
205 92095 : if ((err = mp_mul(&res, &M[bitbuf], &res)) != MP_OKAY) goto LBL_RES;
206 92095 : if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
207 :
208 : /* empty window and reset */
209 87979 : bitcpy = 0;
210 87979 : bitbuf = 0;
211 87979 : mode = 1;
212 : }
213 : }
214 :
215 : /* if bits remain then square/multiply */
216 922 : if ((mode == 2) && (bitcpy > 0)) {
217 : /* square then multiply if the bit is set */
218 2533 : for (x = 0; x < bitcpy; x++) {
219 1655 : if ((err = mp_sqr(&res, &res)) != MP_OKAY) goto LBL_RES;
220 1655 : if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
221 :
222 : /* get next bit of the window */
223 1655 : bitbuf <<= 1;
224 1655 : if ((bitbuf & (1 << winsize)) != 0) {
225 : /* then multiply */
226 1222 : if ((err = mp_mul(&res, &M[1], &res)) != MP_OKAY) goto LBL_RES;
227 1222 : if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
228 : }
229 : }
230 : }
231 :
232 922 : if (redmode == 0) {
233 : /* fixup result if Montgomery reduction is used
234 : * recall that any value in a Montgomery system is
235 : * actually multiplied by R mod n. So we have
236 : * to reduce one more time to cancel out the factor
237 : * of R.
238 : */
239 922 : if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
240 : }
241 :
242 : /* swap res with Y */
243 922 : mp_exch(&res, Y);
244 922 : err = MP_OKAY;
245 922 : LBL_RES:
246 922 : mp_clear(&res);
247 922 : LBL_M:
248 922 : mp_clear(&M[1]);
249 24646 : for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
250 23692 : mp_clear(&M[x]);
251 : }
252 890 : return err;
253 : }
254 : #endif
|
