Mercurial > dropbear

comparison libtommath/bn_mp_exptmod_fast.c @ 1655:f52919ffd3b1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update ltm to 1.1.0 and enable FIPS 186.4 compliant key-generation (#79) * make key-generation compliant to FIPS 186.4 * fix includes in tommath_class.h * update fuzzcorpus instead of error-out * fixup fuzzing make-targets * update Makefile.in * apply necessary patches to ltm sources * clean-up not required ltm files * update to vanilla ltm 1.1.0 this already only contains the required files * remove set/get double
author Steffen Jaeckel <s_jaeckel@gmx.de>
date Mon, 16 Sep 2019 15:50:38 +0200
parents 8bba51a55704
children
comparison
equal deleted inserted replaced
1654:cc0fc5131c5c 1655:f52919ffd3b1
1 #include <tommath_private.h> 1 #include "tommath_private.h"
2 #ifdef BN_MP_EXPTMOD_FAST_C 2 #ifdef BN_MP_EXPTMOD_FAST_C
3 /* LibTomMath, multiple-precision integer library -- Tom St Denis 3 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4 * 4 *
5 * LibTomMath is a library that provides multiple-precision 5 * LibTomMath is a library that provides multiple-precision
6 * integer arithmetic as well as number theoretic functionality. 6 * integer arithmetic as well as number theoretic functionality.
7 * 7 *
8 * The library was designed directly after the MPI library by 8 * The library was designed directly after the MPI library by
9 * Michael Fromberger but has been written from scratch with 9 * Michael Fromberger but has been written from scratch with
10 * additional optimizations in place. 10 * additional optimizations in place.
11 * 11 *
12 * The library is free for all purposes without any express 12 * SPDX-License-Identifier: Unlicense
13 * guarantee it works.
14 *
15 * Tom St Denis, [email protected], http://libtom.org
16 */ 13 */
17 14
18 /* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85 15 /* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
19 * 16 *
20 * Uses a left-to-right k-ary sliding window to compute the modular exponentiation. 17 * Uses a left-to-right k-ary sliding window to compute the modular exponentiation.
22 * 19 *
23 * Uses Montgomery or Diminished Radix reduction [whichever appropriate] 20 * Uses Montgomery or Diminished Radix reduction [whichever appropriate]
24 */ 21 */
25 22
26 #ifdef MP_LOW_MEM 23 #ifdef MP_LOW_MEM
27 #define TAB_SIZE 32 24 # define TAB_SIZE 32
28 #else 25 #else
29 #define TAB_SIZE 256 26 # define TAB_SIZE 256
30 #endif 27 #endif
31 28
32 int mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode) 29 int mp_exptmod_fast(const mp_int *G, const mp_int *X, const mp_int *P, mp_int *Y, int redmode)
33 { 30 {
34 mp_int M[TAB_SIZE], res; 31 mp_int M[TAB_SIZE], res;
35 mp_digit buf, mp; 32 mp_digit buf, mp;
36 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize; 33 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
37 34
38 /* use a pointer to the reduction algorithm. This allows us to use 35 /* use a pointer to the reduction algorithm. This allows us to use
39 * one of many reduction algorithms without modding the guts of 36 * one of many reduction algorithms without modding the guts of
40 * the code with if statements everywhere. 37 * the code with if statements everywhere.
41 */ 38 */
42 int (*redux)(mp_int*,mp_int*,mp_digit); 39 int (*redux)(mp_int *x, const mp_int *n, mp_digit rho);
43 40
44 /* find window size */ 41 /* find window size */
45 x = mp_count_bits (X); 42 x = mp_count_bits(X);
46 if (x <= 7) { 43 if (x <= 7) {
47 winsize = 2; 44 winsize = 2;
48 } else if (x <= 36) { 45 } else if (x <= 36) {
49 winsize = 3; 46 winsize = 3;
50 } else if (x <= 140) { 47 } else if (x <= 140) {
51 winsize = 4; 48 winsize = 4;
52 } else if (x <= 450) { 49 } else if (x <= 450) {
53 winsize = 5; 50 winsize = 5;
54 } else if (x <= 1303) { 51 } else if (x <= 1303) {
55 winsize = 6; 52 winsize = 6;
56 } else if (x <= 3529) { 53 } else if (x <= 3529) {
57 winsize = 7; 54 winsize = 7;
58 } else { 55 } else {
59 winsize = 8; 56 winsize = 8;
60 } 57 }
61 58
62 #ifdef MP_LOW_MEM 59 #ifdef MP_LOW_MEM
63 if (winsize > 5) { 60 if (winsize > 5) {
64 winsize = 5; 61 winsize = 5;
65 } 62 }
66 #endif 63 #endif
67 64
68 /* init M array */ 65 /* init M array */
69 /* init first cell */ 66 /* init first cell */
70 if ((err = mp_init_size(&M[1], P->alloc)) != MP_OKAY) { 67 if ((err = mp_init_size(&M[1], P->alloc)) != MP_OKAY) {
71 return err;
72 }
73
74 /* now init the second half of the array */
75 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
76 if ((err = mp_init_size(&M[x], P->alloc)) != MP_OKAY) {
77 for (y = 1<<(winsize-1); y < x; y++) {
78 mp_clear (&M[y]);
79 }
80 mp_clear(&M[1]);
81 return err; 68 return err;
82 } 69 }
83 } 70
84 71 /* now init the second half of the array */
85 /* determine and setup reduction code */ 72 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
86 if (redmode == 0) { 73 if ((err = mp_init_size(&M[x], P->alloc)) != MP_OKAY) {
87 #ifdef BN_MP_MONTGOMERY_SETUP_C 74 for (y = 1<<(winsize-1); y < x; y++) {
88 /* now setup montgomery */ 75 mp_clear(&M[y]);
89 if ((err = mp_montgomery_setup (P, &mp)) != MP_OKAY) { 76 }
90 goto LBL_M; 77 mp_clear(&M[1]);
91 } 78 return err;
92 #else 79 }
93 err = MP_VAL; 80 }
94 goto LBL_M; 81
95 #endif 82 /* determine and setup reduction code */
96 83 if (redmode == 0) {
97 /* automatically pick the comba one if available (saves quite a few calls/ifs) */ 84 #ifdef BN_MP_MONTGOMERY_SETUP_C
85 /* now setup montgomery */
86 if ((err = mp_montgomery_setup(P, &mp)) != MP_OKAY) {
87 goto LBL_M;
88 }
89 #else
90 err = MP_VAL;
91 goto LBL_M;
92 #endif
93
94 /* automatically pick the comba one if available (saves quite a few calls/ifs) */
98 #ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C 95 #ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C
99 if ((((P->used * 2) + 1) < MP_WARRAY) && 96 if ((((P->used * 2) + 1) < (int)MP_WARRAY) &&
100 (P->used < (1 << ((CHAR_BIT * sizeof(mp_word)) - (2 * DIGIT_BIT))))) { 97 (P->used < (1 << ((CHAR_BIT * sizeof(mp_word)) - (2 * DIGIT_BIT))))) {
101 redux = fast_mp_montgomery_reduce; 98 redux = fast_mp_montgomery_reduce;
102 } else 99 } else
103 #endif 100 #endif
104 { 101 {
105 #ifdef BN_MP_MONTGOMERY_REDUCE_C 102 #ifdef BN_MP_MONTGOMERY_REDUCE_C
106 /* use slower baseline Montgomery method */ 103 /* use slower baseline Montgomery method */
107 redux = mp_montgomery_reduce; 104 redux = mp_montgomery_reduce;
108 #else 105 #else
109 err = MP_VAL; 106 err = MP_VAL;
110 goto LBL_M; 107 goto LBL_M;
111 #endif 108 #endif
112 } 109 }
113 } else if (redmode == 1) { 110 } else if (redmode == 1) {
114 #if defined(BN_MP_DR_SETUP_C) && defined(BN_MP_DR_REDUCE_C) 111 #if defined(BN_MP_DR_SETUP_C) && defined(BN_MP_DR_REDUCE_C)
115 /* setup DR reduction for moduli of the form B**k - b */ 112 /* setup DR reduction for moduli of the form B**k - b */
116 mp_dr_setup(P, &mp); 113 mp_dr_setup(P, &mp);
117 redux = mp_dr_reduce; 114 redux = mp_dr_reduce;
118 #else 115 #else
119 err = MP_VAL; 116 err = MP_VAL;
120 goto LBL_M; 117 goto LBL_M;
121 #endif 118 #endif
122 } else { 119 } else {
123 #if defined(BN_MP_REDUCE_2K_SETUP_C) && defined(BN_MP_REDUCE_2K_C) 120 #if defined(BN_MP_REDUCE_2K_SETUP_C) && defined(BN_MP_REDUCE_2K_C)
124 /* setup DR reduction for moduli of the form 2**k - b */ 121 /* setup DR reduction for moduli of the form 2**k - b */
125 if ((err = mp_reduce_2k_setup(P, &mp)) != MP_OKAY) { 122 if ((err = mp_reduce_2k_setup(P, &mp)) != MP_OKAY) {
126 goto LBL_M; 123 goto LBL_M;
127 } 124 }
128 redux = mp_reduce_2k; 125 redux = mp_reduce_2k;
129 #else 126 #else
130 err = MP_VAL; 127 err = MP_VAL;
131 goto LBL_M; 128 goto LBL_M;
132 #endif 129 #endif
133 } 130 }
134 131
135 /* setup result */ 132 /* setup result */
136 if ((err = mp_init_size (&res, P->alloc)) != MP_OKAY) { 133 if ((err = mp_init_size(&res, P->alloc)) != MP_OKAY) {
137 goto LBL_M; 134 goto LBL_M;
138 } 135 }
139 136
140 /* create M table 137 /* create M table
141 * 138 *
142 139
143 * 140 *
144 * The first half of the table is not computed though accept for M[0] and M[1] 141 * The first half of the table is not computed though accept for M[0] and M[1]
145 */ 142 */
146 143
147 if (redmode == 0) { 144 if (redmode == 0) {
148 #ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C 145 #ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C
149 /* now we need R mod m */ 146 /* now we need R mod m */
150 if ((err = mp_montgomery_calc_normalization (&res, P)) != MP_OKAY) { 147 if ((err = mp_montgomery_calc_normalization(&res, P)) != MP_OKAY) {
151 goto LBL_RES; 148 goto LBL_RES;
152 } 149 }
153 150
154 /* now set M[1] to G * R mod m */ 151 /* now set M[1] to G * R mod m */
155 if ((err = mp_mulmod (G, &res, P, &M[1])) != MP_OKAY) { 152 if ((err = mp_mulmod(G, &res, P, &M[1])) != MP_OKAY) {
156 goto LBL_RES; 153 goto LBL_RES;
157 } 154 }
158 #else 155 #else
159 err = MP_VAL; 156 err = MP_VAL;
160 goto LBL_RES;
161 #endif
162 } else {
163 mp_set(&res, 1);
164 if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {
165 goto LBL_RES;
166 }
167 }
168
169 /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
170 if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
171 goto LBL_RES;
172 }
173
174 for (x = 0; x < (winsize - 1); x++) {
175 if ((err = mp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) {
176 goto LBL_RES; 157 goto LBL_RES;
177 } 158 #endif
178 if ((err = redux (&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) { 159 } else {
160 mp_set(&res, 1uL);
161 if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {
162 goto LBL_RES;
163 }
164 }
165
166 /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
167 if ((err = mp_copy(&M[1], &M[(size_t)1 << (winsize - 1)])) != MP_OKAY) {
179 goto LBL_RES; 168 goto LBL_RES;
180 } 169 }
181 } 170
182 171 for (x = 0; x < (winsize - 1); x++) {
183 /* create upper table */ 172 if ((err = mp_sqr(&M[(size_t)1 << (winsize - 1)], &M[(size_t)1 << (winsize - 1)])) != MP_OKAY) {
184 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) { 173 goto LBL_RES;
185 if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) { 174 }
186 goto LBL_RES; 175 if ((err = redux(&M[(size_t)1 << (winsize - 1)], P, mp)) != MP_OKAY) {
187 } 176 goto LBL_RES;
188 if ((err = redux (&M[x], P, mp)) != MP_OKAY) { 177 }
189 goto LBL_RES; 178 }
190 } 179
191 } 180 /* create upper table */
192 181 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
193 /* set initial mode and bit cnt */ 182 if ((err = mp_mul(&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
194 mode = 0; 183 goto LBL_RES;
195 bitcnt = 1; 184 }
196 buf = 0; 185 if ((err = redux(&M[x], P, mp)) != MP_OKAY) {
197 digidx = X->used - 1; 186 goto LBL_RES;
198 bitcpy = 0; 187 }
199 bitbuf = 0; 188 }
200 189
201 for (;;) { 190 /* set initial mode and bit cnt */
202 /* grab next digit as required */ 191 mode = 0;
203 if (--bitcnt == 0) { 192 bitcnt = 1;
204 /* if digidx == -1 we are out of digits so break */ 193 buf = 0;
205 if (digidx == -1) { 194 digidx = X->used - 1;
206 break; 195 bitcpy = 0;
207 } 196 bitbuf = 0;
208 /* read next digit and reset bitcnt */ 197
209 buf = X->dp[digidx--]; 198 for (;;) {
210 bitcnt = (int)DIGIT_BIT; 199 /* grab next digit as required */
211 } 200 if (--bitcnt == 0) {
212 201 /* if digidx == -1 we are out of digits so break */
213 /* grab the next msb from the exponent */ 202 if (digidx == -1) {
214 y = (mp_digit)(buf >> (DIGIT_BIT - 1)) & 1; 203 break;
215 buf <<= (mp_digit)1; 204 }
216 205 /* read next digit and reset bitcnt */
217 /* if the bit is zero and mode == 0 then we ignore it 206 buf = X->dp[digidx--];
218 * These represent the leading zero bits before the first 1 bit 207 bitcnt = (int)DIGIT_BIT;
219 * in the exponent. Technically this opt is not required but it 208 }
220 * does lower the # of trivial squaring/reductions used 209
221 */ 210 /* grab the next msb from the exponent */
222 if ((mode == 0) && (y == 0)) { 211 y = (mp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
223 continue; 212 buf <<= (mp_digit)1;
224 } 213
225 214 /* if the bit is zero and mode == 0 then we ignore it
226 /* if the bit is zero and mode == 1 then we square */ 215 * These represent the leading zero bits before the first 1 bit
227 if ((mode == 1) && (y == 0)) { 216 * in the exponent. Technically this opt is not required but it
228 if ((err = mp_sqr (&res, &res)) != MP_OKAY) { 217 * does lower the # of trivial squaring/reductions used
229 goto LBL_RES; 218 */
230 } 219 if ((mode == 0) && (y == 0)) {
231 if ((err = redux (&res, P, mp)) != MP_OKAY) { 220 continue;
232 goto LBL_RES; 221 }
233 } 222
234 continue; 223 /* if the bit is zero and mode == 1 then we square */
235 } 224 if ((mode == 1) && (y == 0)) {
236 225 if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
237 /* else we add it to the window */ 226 goto LBL_RES;
238 bitbuf |= (y << (winsize - ++bitcpy)); 227 }
239 mode = 2; 228 if ((err = redux(&res, P, mp)) != MP_OKAY) {
240 229 goto LBL_RES;
241 if (bitcpy == winsize) { 230 }
242 /* ok window is filled so square as required and multiply */ 231 continue;
243 /* square first */ 232 }
244 for (x = 0; x < winsize; x++) { 233
245 if ((err = mp_sqr (&res, &res)) != MP_OKAY) { 234 /* else we add it to the window */
246 goto LBL_RES; 235 bitbuf |= (y << (winsize - ++bitcpy));
247 } 236 mode = 2;
248 if ((err = redux (&res, P, mp)) != MP_OKAY) { 237
249 goto LBL_RES; 238 if (bitcpy == winsize) {
250 } 239 /* ok window is filled so square as required and multiply */
251 } 240 /* square first */
252 241 for (x = 0; x < winsize; x++) {
253 /* then multiply */ 242 if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
254 if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) { 243 goto LBL_RES;
255 goto LBL_RES; 244 }
256 } 245 if ((err = redux(&res, P, mp)) != MP_OKAY) {
257 if ((err = redux (&res, P, mp)) != MP_OKAY) { 246 goto LBL_RES;
258 goto LBL_RES; 247 }
259 } 248 }
260 249
261 /* empty window and reset */ 250 /* then multiply */
262 bitcpy = 0; 251 if ((err = mp_mul(&res, &M[bitbuf], &res)) != MP_OKAY) {
263 bitbuf = 0; 252 goto LBL_RES;
264 mode = 1; 253 }
265 } 254 if ((err = redux(&res, P, mp)) != MP_OKAY) {
266 } 255 goto LBL_RES;
267 256 }
268 /* if bits remain then square/multiply */ 257
269 if ((mode == 2) && (bitcpy > 0)) { 258 /* empty window and reset */
270 /* square then multiply if the bit is set */ 259 bitcpy = 0;
271 for (x = 0; x < bitcpy; x++) { 260 bitbuf = 0;
272 if ((err = mp_sqr (&res, &res)) != MP_OKAY) { 261 mode = 1;
273 goto LBL_RES; 262 }
274 } 263 }
275 if ((err = redux (&res, P, mp)) != MP_OKAY) { 264
276 goto LBL_RES; 265 /* if bits remain then square/multiply */
277 } 266 if ((mode == 2) && (bitcpy > 0)) {
278 267 /* square then multiply if the bit is set */
279 /* get next bit of the window */ 268 for (x = 0; x < bitcpy; x++) {
280 bitbuf <<= 1; 269 if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
281 if ((bitbuf & (1 << winsize)) != 0) { 270 goto LBL_RES;
282 /* then multiply */ 271 }
283 if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) { 272 if ((err = redux(&res, P, mp)) != MP_OKAY) {
284 goto LBL_RES; 273 goto LBL_RES;
285 } 274 }
286 if ((err = redux (&res, P, mp)) != MP_OKAY) { 275
287 goto LBL_RES; 276 /* get next bit of the window */
288 } 277 bitbuf <<= 1;
289 } 278 if ((bitbuf & (1 << winsize)) != 0) {
290 } 279 /* then multiply */
291 } 280 if ((err = mp_mul(&res, &M[1], &res)) != MP_OKAY) {
292 281 goto LBL_RES;
293 if (redmode == 0) { 282 }
294 /* fixup result if Montgomery reduction is used 283 if ((err = redux(&res, P, mp)) != MP_OKAY) {
295 * recall that any value in a Montgomery system is 284 goto LBL_RES;
296 * actually multiplied by R mod n. So we have 285 }
297 * to reduce one more time to cancel out the factor 286 }
298 * of R. 287 }
299 */ 288 }
300 if ((err = redux(&res, P, mp)) != MP_OKAY) { 289
301 goto LBL_RES; 290 if (redmode == 0) {
302 } 291 /* fixup result if Montgomery reduction is used
303 } 292 * recall that any value in a Montgomery system is
304 293 * actually multiplied by R mod n. So we have
305 /* swap res with Y */ 294 * to reduce one more time to cancel out the factor
306 mp_exch (&res, Y); 295 * of R.
307 err = MP_OKAY; 296 */
308 LBL_RES:mp_clear (&res); 297 if ((err = redux(&res, P, mp)) != MP_OKAY) {
298 goto LBL_RES;
299 }
300 }
301
302 /* swap res with Y */
303 mp_exch(&res, Y);
304 err = MP_OKAY;
305 LBL_RES:
306 mp_clear(&res);
309 LBL_M: 307 LBL_M:
310 mp_clear(&M[1]); 308 mp_clear(&M[1]);
311 for (x = 1<<(winsize-1); x < (1 << winsize); x++) { 309 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
312 mp_clear (&M[x]); 310 mp_clear(&M[x]);
313 } 311 }
314 return err; 312 return err;
315 } 313 }
316 #endif 314 #endif
317 315
318 316
319 /* ref: $Format:%D$ */ 317 /* ref: HEAD -> master, tag: v1.1.0 */
320 /* git commit: $Format:%H$ */ 318 /* git commit: 08549ad6bc8b0cede0b357a9c341c5c6473a9c55 */
321 /* commit time: $Format:%ai$ */ 319 /* commit time: 2019-01-28 20:32:32 +0100 */