1 #include "tommath_private.h"
2 #ifdef BN_S_MP_EXPTMOD_C
3 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4 *
5 * LibTomMath is a library that provides multiple-precision
6 * integer arithmetic as well as number theoretic functionality.
7 *
8 * The library was designed directly after the MPI library by
9 * Michael Fromberger but has been written from scratch with
10 * additional optimizations in place.
11 *
12 * SPDX-License-Identifier: Unlicense
13 */
14
15 #ifdef MP_LOW_MEM
16 # define TAB_SIZE 32
17 #else
18 # define TAB_SIZE 256
19 #endif
20
s_mp_exptmod(const mp_int * G,const mp_int * X,const mp_int * P,mp_int * Y,int redmode)21 int s_mp_exptmod(const mp_int *G, const mp_int *X, const mp_int *P, mp_int *Y, int redmode)
22 {
23 mp_int M[TAB_SIZE], res, mu;
24 mp_digit buf;
25 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
26 int (*redux)(mp_int *x, const mp_int *m, const mp_int *mu);
27
28 /* find window size */
29 x = mp_count_bits(X);
30 if (x <= 7) {
31 winsize = 2;
32 } else if (x <= 36) {
33 winsize = 3;
34 } else if (x <= 140) {
35 winsize = 4;
36 } else if (x <= 450) {
37 winsize = 5;
38 } else if (x <= 1303) {
39 winsize = 6;
40 } else if (x <= 3529) {
41 winsize = 7;
42 } else {
43 winsize = 8;
44 }
45
46 #ifdef MP_LOW_MEM
47 if (winsize > 5) {
48 winsize = 5;
49 }
50 #endif
51
52 /* init M array */
53 /* init first cell */
54 if ((err = mp_init(&M[1])) != MP_OKAY) {
55 return err;
56 }
57
58 /* now init the second half of the array */
59 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
60 if ((err = mp_init(&M[x])) != MP_OKAY) {
61 for (y = 1<<(winsize-1); y < x; y++) {
62 mp_clear(&M[y]);
63 }
64 mp_clear(&M[1]);
65 return err;
66 }
67 }
68
69 /* create mu, used for Barrett reduction */
70 if ((err = mp_init(&mu)) != MP_OKAY) {
71 goto LBL_M;
72 }
73
74 if (redmode == 0) {
75 if ((err = mp_reduce_setup(&mu, P)) != MP_OKAY) {
76 goto LBL_MU;
77 }
78 redux = mp_reduce;
79 } else {
80 if ((err = mp_reduce_2k_setup_l(P, &mu)) != MP_OKAY) {
81 goto LBL_MU;
82 }
83 redux = mp_reduce_2k_l;
84 }
85
86 /* create M table
87 *
88 * The M table contains powers of the base,
89 * e.g. M[x] = G**x mod P
90 *
91 * The first half of the table is not
92 * computed though accept for M[0] and M[1]
93 */
94 if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {
95 goto LBL_MU;
96 }
97
98 /* compute the value at M[1<<(winsize-1)] by squaring
99 * M[1] (winsize-1) times
100 */
101 if ((err = mp_copy(&M[1], &M[(size_t)1 << (winsize - 1)])) != MP_OKAY) {
102 goto LBL_MU;
103 }
104
105 for (x = 0; x < (winsize - 1); x++) {
106 /* square it */
107 if ((err = mp_sqr(&M[(size_t)1 << (winsize - 1)],
108 &M[(size_t)1 << (winsize - 1)])) != MP_OKAY) {
109 goto LBL_MU;
110 }
111
112 /* reduce modulo P */
113 if ((err = redux(&M[(size_t)1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
114 goto LBL_MU;
115 }
116 }
117
118 /* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
119 * for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
120 */
121 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
122 if ((err = mp_mul(&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
123 goto LBL_MU;
124 }
125 if ((err = redux(&M[x], P, &mu)) != MP_OKAY) {
126 goto LBL_MU;
127 }
128 }
129
130 /* setup result */
131 if ((err = mp_init(&res)) != MP_OKAY) {
132 goto LBL_MU;
133 }
134 mp_set(&res, 1uL);
135
136 /* set initial mode and bit cnt */
137 mode = 0;
138 bitcnt = 1;
139 buf = 0;
140 digidx = X->used - 1;
141 bitcpy = 0;
142 bitbuf = 0;
143
144 for (;;) {
145 /* grab next digit as required */
146 if (--bitcnt == 0) {
147 /* if digidx == -1 we are out of digits */
148 if (digidx == -1) {
149 break;
150 }
151 /* read next digit and reset the bitcnt */
152 buf = X->dp[digidx--];
153 bitcnt = (int)DIGIT_BIT;
154 }
155
156 /* grab the next msb from the exponent */
157 y = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1;
158 buf <<= (mp_digit)1;
159
160 /* if the bit is zero and mode == 0 then we ignore it
161 * These represent the leading zero bits before the first 1 bit
162 * in the exponent. Technically this opt is not required but it
163 * does lower the # of trivial squaring/reductions used
164 */
165 if ((mode == 0) && (y == 0)) {
166 continue;
167 }
168
169 /* if the bit is zero and mode == 1 then we square */
170 if ((mode == 1) && (y == 0)) {
171 if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
172 goto LBL_RES;
173 }
174 if ((err = redux(&res, P, &mu)) != MP_OKAY) {
175 goto LBL_RES;
176 }
177 continue;
178 }
179
180 /* else we add it to the window */
181 bitbuf |= (y << (winsize - ++bitcpy));
182 mode = 2;
183
184 if (bitcpy == winsize) {
185 /* ok window is filled so square as required and multiply */
186 /* square first */
187 for (x = 0; x < winsize; x++) {
188 if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
189 goto LBL_RES;
190 }
191 if ((err = redux(&res, P, &mu)) != MP_OKAY) {
192 goto LBL_RES;
193 }
194 }
195
196 /* then multiply */
197 if ((err = mp_mul(&res, &M[bitbuf], &res)) != MP_OKAY) {
198 goto LBL_RES;
199 }
200 if ((err = redux(&res, P, &mu)) != MP_OKAY) {
201 goto LBL_RES;
202 }
203
204 /* empty window and reset */
205 bitcpy = 0;
206 bitbuf = 0;
207 mode = 1;
208 }
209 }
210
211 /* if bits remain then square/multiply */
212 if ((mode == 2) && (bitcpy > 0)) {
213 /* square then multiply if the bit is set */
214 for (x = 0; x < bitcpy; x++) {
215 if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
216 goto LBL_RES;
217 }
218 if ((err = redux(&res, P, &mu)) != MP_OKAY) {
219 goto LBL_RES;
220 }
221
222 bitbuf <<= 1;
223 if ((bitbuf & (1 << winsize)) != 0) {
224 /* then multiply */
225 if ((err = mp_mul(&res, &M[1], &res)) != MP_OKAY) {
226 goto LBL_RES;
227 }
228 if ((err = redux(&res, P, &mu)) != MP_OKAY) {
229 goto LBL_RES;
230 }
231 }
232 }
233 }
234
235 mp_exch(&res, Y);
236 err = MP_OKAY;
237 LBL_RES:
238 mp_clear(&res);
239 LBL_MU:
240 mp_clear(&mu);
241 LBL_M:
242 mp_clear(&M[1]);
243 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
244 mp_clear(&M[x]);
245 }
246 return err;
247 }
248 #endif
249
250 /* ref: $Format:%D$ */
251 /* git commit: $Format:%H$ */
252 /* commit time: $Format:%ai$ */
253