Mercurial > dropbear

comparison libtommath/bn_s_mp_exptmod_fast.c @ 1692:1051e4eea25a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update LibTomMath to 1.2.0 (#84) * update C files * update other files * update headers * update makefiles * remove mp_set/get_double() * use ltm 1.2.0 API * update ltm_desc * use bundled tommath if system-tommath is too old * XMALLOC etc. were changed to MP_MALLOC etc.
author Steffen Jaeckel <s@jaeckel.eu>
date Tue, 26 May 2020 17:36:47 +0200
parents
children
comparison
equal deleted inserted replaced
1691:2d3745d58843 1692:1051e4eea25a
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 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 mp_int M[TAB_SIZE], res;
25 mp_digit buf, mp;
26 int bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
27 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 mp_err(*redux)(mp_int *x, const mp_int *n, mp_digit rho);
34
35 /* find window size */
36 x = mp_count_bits(X);
37 if (x <= 7) {
38 winsize = 2;
39 } else if (x <= 36) {
40 winsize = 3;
41 } else if (x <= 140) {
42 winsize = 4;
43 } else if (x <= 450) {
44 winsize = 5;
45 } else if (x <= 1303) {
46 winsize = 6;
47 } else if (x <= 3529) {
48 winsize = 7;
49 } else {
50 winsize = 8;
51 }
52
53 winsize = MAX_WINSIZE ? MP_MIN(MAX_WINSIZE, winsize) : winsize;
54
55 /* init M array */
56 /* init first cell */
57 if ((err = mp_init_size(&M[1], P->alloc)) != MP_OKAY) {
58 return err;
59 }
60
61 /* now init the second half of the array */
62 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
63 if ((err = mp_init_size(&M[x], P->alloc)) != MP_OKAY) {
64 for (y = 1<<(winsize-1); y < x; y++) {
65 mp_clear(&M[y]);
66 }
67 mp_clear(&M[1]);
68 return err;
69 }
70 }
71
72 /* determine and setup reduction code */
73 if (redmode == 0) {
74 if (MP_HAS(MP_MONTGOMERY_SETUP)) {
75 /* now setup montgomery */
76 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 if (MP_HAS(S_MP_MONTGOMERY_REDUCE_FAST) &&
84 (((P->used * 2) + 1) < MP_WARRAY) &&
85 (P->used < MP_MAXFAST)) {
86 redux = s_mp_montgomery_reduce_fast;
87 } else if (MP_HAS(MP_MONTGOMERY_REDUCE)) {
88 /* use slower baseline Montgomery method */
89 redux = mp_montgomery_reduce;
90 } else {
91 err = MP_VAL;
92 goto LBL_M;
93 }
94 } else if (redmode == 1) {
95 if (MP_HAS(MP_DR_SETUP) && MP_HAS(MP_DR_REDUCE)) {
96 /* setup DR reduction for moduli of the form B**k - b */
97 mp_dr_setup(P, &mp);
98 redux = mp_dr_reduce;
99 } else {
100 err = MP_VAL;
101 goto LBL_M;
102 }
103 } 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 if ((err = mp_reduce_2k_setup(P, &mp)) != MP_OKAY) goto LBL_M;
106 redux = mp_reduce_2k;
107 } else {
108 err = MP_VAL;
109 goto LBL_M;
110 }
111
112 /* setup result */
113 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 if (redmode == 0) {
123 if (MP_HAS(MP_MONTGOMERY_CALC_NORMALIZATION)) {
124 /* now we need R mod m */
125 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 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 mp_set(&res, 1uL);
135 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 if ((err = mp_copy(&M[1], &M[(size_t)1 << (winsize - 1)])) != MP_OKAY) goto LBL_RES;
140
141 for (x = 0; x < (winsize - 1); x++) {
142 if ((err = mp_sqr(&M[(size_t)1 << (winsize - 1)], &M[(size_t)1 << (winsize - 1)])) != MP_OKAY) goto LBL_RES;
143 if ((err = redux(&M[(size_t)1 << (winsize - 1)], P, mp)) != MP_OKAY) goto LBL_RES;
144 }
145
146 /* create upper table */
147 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
148 if ((err = mp_mul(&M[x - 1], &M[1], &M[x])) != MP_OKAY) goto LBL_RES;
149 if ((err = redux(&M[x], P, mp)) != MP_OKAY) goto LBL_RES;
150 }
151
152 /* set initial mode and bit cnt */
153 mode = 0;
154 bitcnt = 1;
155 buf = 0;
156 digidx = X->used - 1;
157 bitcpy = 0;
158 bitbuf = 0;
159
160 for (;;) {
161 /* grab next digit as required */
162 if (--bitcnt == 0) {
163 /* if digidx == -1 we are out of digits so break */
164 if (digidx == -1) {
165 break;
166 }
167 /* read next digit and reset bitcnt */
168 buf = X->dp[digidx--];
169 bitcnt = (int)MP_DIGIT_BIT;
170 }
171
172 /* grab the next msb from the exponent */
173 y = (mp_digit)(buf >> (MP_DIGIT_BIT - 1)) & 1uL;
174 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 if ((mode == 0) && (y == 0)) {
182 continue;
183 }
184
185 /* if the bit is zero and mode == 1 then we square */
186 if ((mode == 1) && (y == 0)) {
187 if ((err = mp_sqr(&res, &res)) != MP_OKAY) goto LBL_RES;
188 if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
189 continue;
190 }
191
192 /* else we add it to the window */
193 bitbuf |= (y << (winsize - ++bitcpy));
194 mode = 2;
195
196 if (bitcpy == winsize) {
197 /* ok window is filled so square as required and multiply */
198 /* square first */
199 for (x = 0; x < winsize; x++) {
200 if ((err = mp_sqr(&res, &res)) != MP_OKAY) goto LBL_RES;
201 if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
202 }
203
204 /* then multiply */
205 if ((err = mp_mul(&res, &M[bitbuf], &res)) != MP_OKAY) goto LBL_RES;
206 if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
207
208 /* empty window and reset */
209 bitcpy = 0;
210 bitbuf = 0;
211 mode = 1;
212 }
213 }
214
215 /* if bits remain then square/multiply */
216 if ((mode == 2) && (bitcpy > 0)) {
217 /* square then multiply if the bit is set */
218 for (x = 0; x < bitcpy; x++) {
219 if ((err = mp_sqr(&res, &res)) != MP_OKAY) goto LBL_RES;
220 if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
221
222 /* get next bit of the window */
223 bitbuf <<= 1;
224 if ((bitbuf & (1 << winsize)) != 0) {
225 /* then multiply */
226 if ((err = mp_mul(&res, &M[1], &res)) != MP_OKAY) goto LBL_RES;
227 if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
228 }
229 }
230 }
231
232 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 if ((err = redux(&res, P, mp)) != MP_OKAY) goto LBL_RES;
240 }
241
242 /* swap res with Y */
243 mp_exch(&res, Y);
244 err = MP_OKAY;
245 LBL_RES:
246 mp_clear(&res);
247 LBL_M:
248 mp_clear(&M[1]);
249 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
250 mp_clear(&M[x]);
251 }
252 return err;
253 }
254 #endif