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