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 893 : 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 893 : x = mp_count_bits(X);
37 893 : if (x <= 7) {
38 0 : winsize = 2;
39 893 : } else if (x <= 36) {
40 422 : winsize = 3;
41 455 : } else if (x <= 140) {
42 0 : winsize = 4;
43 455 : } else if (x <= 450) {
44 0 : winsize = 5;
45 455 : } else if (x <= 1303) {
46 243 : winsize = 6;
47 212 : } else if (x <= 3529) {
48 196 : winsize = 7;
49 : } else {
50 0 : winsize = 8;
51 : }
52 :
53 893 : winsize = MAX_WINSIZE ? MP_MIN(MAX_WINSIZE, winsize) : winsize;
54 :
55 : /* init M array */
56 : /* init first cell */
57 893 : 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 23989 : for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
63 23096 : 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 893 : if (redmode == 0) {
74 32 : if (MP_HAS(MP_MONTGOMERY_SETUP)) {
75 : /* now setup montgomery */
76 893 : 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 893 : if (MP_HAS(S_MP_MONTGOMERY_REDUCE_FAST) &&
84 893 : (((P->used * 2) + 1) < MP_WARRAY) &&
85 861 : (P->used < MP_MAXFAST)) {
86 861 : 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 893 : 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 893 : if (redmode == 0) {
123 32 : if (MP_HAS(MP_MONTGOMERY_CALC_NORMALIZATION)) {
124 : /* now we need R mod m */
125 893 : 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 893 : 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 893 : if ((err = mp_copy(&M[1], &M[(size_t)1 << (winsize - 1)])) != MP_OKAY) goto LBL_RES;
140 :
141 4256 : for (x = 0; x < (winsize - 1); x++) {
142 3363 : if ((err = mp_sqr(&M[(size_t)1 << (winsize - 1)], &M[(size_t)1 << (winsize - 1)])) != MP_OKAY) goto LBL_RES;
143 3363 : if ((err = redux(&M[(size_t)1 << (winsize - 1)], P, mp)) != MP_OKAY) goto LBL_RES;
144 : }
145 :
146 : /* create upper table */
147 23096 : for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
148 22203 : if ((err = mp_mul(&M[x - 1], &M[1], &M[x])) != MP_OKAY) goto LBL_RES;
149 22203 : if ((err = redux(&M[x], P, mp)) != MP_OKAY) goto LBL_RES;
150 : }
151 :
152 : /* set initial mode and bit cnt */
153 893 : mode = 0;
154 893 : bitcnt = 1;
155 893 : buf = 0;
156 893 : digidx = X->used - 1;
157 893 : bitcpy = 0;
158 893 : bitbuf = 0;
159 :
160 34592 : for (;;) {
161 : /* grab next digit as required */
162 734813 : if (--bitcnt == 0) {
163 : /* if digidx == -1 we are out of digits so break */
164 13125 : if (digidx == -1) {
165 861 : break;
166 : }
167 : /* read next digit and reset bitcnt */
168 12232 : buf = X->dp[digidx--];
169 12232 : bitcnt = (int)MP_DIGIT_BIT;
170 : }
171 :
172 : /* grab the next msb from the exponent */
173 733920 : y = (mp_digit)(buf >> (MP_DIGIT_BIT - 1)) & 1uL;
174 733920 : 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 733920 : if ((mode == 0) && (y == 0)) {
182 43981 : continue;
183 : }
184 :
185 : /* if the bit is zero and mode == 1 then we square */
186 689939 : if ((mode == 1) && (y == 0)) {
187 94967 : if ((err = mp_sqr(&res, &res)) != MP_OKAY) goto LBL_RES;
188 94967 : if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
189 94967 : continue;
190 : }
191 :
192 : /* else we add it to the window */
193 594972 : bitbuf |= (y << (winsize - ++bitcpy));
194 594972 : mode = 2;
195 :
196 594972 : if (bitcpy == winsize) {
197 : /* ok window is filled so square as required and multiply */
198 : /* square first */
199 683506 : for (x = 0; x < winsize; x++) {
200 593415 : if ((err = mp_sqr(&res, &res)) != MP_OKAY) goto LBL_RES;
201 593415 : if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
202 : }
203 :
204 : /* then multiply */
205 90091 : if ((err = mp_mul(&res, &M[bitbuf], &res)) != MP_OKAY) goto LBL_RES;
206 90091 : if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
207 :
208 : /* empty window and reset */
209 85975 : bitcpy = 0;
210 85975 : bitbuf = 0;
211 85975 : mode = 1;
212 : }
213 : }
214 :
215 : /* if bits remain then square/multiply */
216 893 : if ((mode == 2) && (bitcpy > 0)) {
217 : /* square then multiply if the bit is set */
218 2407 : for (x = 0; x < bitcpy; x++) {
219 1557 : if ((err = mp_sqr(&res, &res)) != MP_OKAY) goto LBL_RES;
220 1557 : if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
221 :
222 : /* get next bit of the window */
223 1557 : bitbuf <<= 1;
224 1557 : if ((bitbuf & (1 << winsize)) != 0) {
225 : /* then multiply */
226 1167 : if ((err = mp_mul(&res, &M[1], &res)) != MP_OKAY) goto LBL_RES;
227 1167 : if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
228 : }
229 : }
230 : }
231 :
232 893 : 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 893 : if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
240 : }
241 :
242 : /* swap res with Y */
243 893 : mp_exch(&res, Y);
244 893 : err = MP_OKAY;
245 893 : LBL_RES:
246 893 : mp_clear(&res);
247 893 : LBL_M:
248 893 : mp_clear(&M[1]);
249 24021 : for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
250 23096 : mp_clear(&M[x]);
251 : }
252 861 : return err;
253 : }
254 : #endif
|
