Mercurial > dropbear

comparison bn_s_mp_exptmod.c @ 282:91fbc376f010 libtommath-orig libtommath-0.35

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import of libtommath 0.35 From ltm-0.35.tar.bz2 SHA1 of 3f193dbae9351e92d02530994fa18236f7fde01c
author Matt Johnston <matt@ucc.asn.au>
date Wed, 08 Mar 2006 13:16:18 +0000
parents
children 97db060d0ef5
comparison
equal deleted inserted replaced
-1:000000000000 282:91fbc376f010
1 #include <tommath.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 * The library is free for all purposes without any express
13 * guarantee it works.
14 *
15 * Tom St Denis, [email protected], http://math.libtomcrypt.org
16 */
17
18 #ifdef MP_LOW_MEM
19 #define TAB_SIZE 32
20 #else
21 #define TAB_SIZE 256
22 #endif
23
24 int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
25 {
26 mp_int M[TAB_SIZE], res, mu;
27 mp_digit buf;
28 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
29 int (*redux)(mp_int*,mp_int*,mp_int*);
30
31 /* find window size */
32 x = mp_count_bits (X);
33 if (x <= 7) {
34 winsize = 2;
35 } else if (x <= 36) {
36 winsize = 3;
37 } else if (x <= 140) {
38 winsize = 4;
39 } else if (x <= 450) {
40 winsize = 5;
41 } else if (x <= 1303) {
42 winsize = 6;
43 } else if (x <= 3529) {
44 winsize = 7;
45 } else {
46 winsize = 8;
47 }
48
49 #ifdef MP_LOW_MEM
50 if (winsize > 5) {
51 winsize = 5;
52 }
53 #endif
54
55 /* init M array */
56 /* init first cell */
57 if ((err = mp_init(&M[1])) != 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(&M[x])) != 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 /* create mu, used for Barrett reduction */
73 if ((err = mp_init (&mu)) != MP_OKAY) {
74 goto LBL_M;
75 }
76
77 if (redmode == 0) {
78 if ((err = mp_reduce_setup (&mu, P)) != MP_OKAY) {
79 goto LBL_MU;
80 }
81 redux = mp_reduce;
82 } else {
83 if ((err = mp_reduce_2k_setup_l (P, &mu)) != MP_OKAY) {
84 goto LBL_MU;
85 }
86 redux = mp_reduce_2k_l;
87 }
88
89 /* create M table
90 *
91 * The M table contains powers of the base,
92 * e.g. M[x] = G**x mod P
93 *
94 * The first half of the table is not
95 * computed though accept for M[0] and M[1]
96 */
97 if ((err = mp_mod (G, P, &M[1])) != MP_OKAY) {
98 goto LBL_MU;
99 }
100
101 /* compute the value at M[1<<(winsize-1)] by squaring
102 * M[1] (winsize-1) times
103 */
104 if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
105 goto LBL_MU;
106 }
107
108 for (x = 0; x < (winsize - 1); x++) {
109 /* square it */
110 if ((err = mp_sqr (&M[1 << (winsize - 1)],
111 &M[1 << (winsize - 1)])) != MP_OKAY) {
112 goto LBL_MU;
113 }
114
115 /* reduce modulo P */
116 if ((err = redux (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
117 goto LBL_MU;
118 }
119 }
120
121 /* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
122 * for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
123 */
124 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
125 if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
126 goto LBL_MU;
127 }
128 if ((err = redux (&M[x], P, &mu)) != MP_OKAY) {
129 goto LBL_MU;
130 }
131 }
132
133 /* setup result */
134 if ((err = mp_init (&res)) != MP_OKAY) {
135 goto LBL_MU;
136 }
137 mp_set (&res, 1);
138
139 /* set initial mode and bit cnt */
140 mode = 0;
141 bitcnt = 1;
142 buf = 0;
143 digidx = X->used - 1;
144 bitcpy = 0;
145 bitbuf = 0;
146
147 for (;;) {
148 /* grab next digit as required */
149 if (--bitcnt == 0) {
150 /* if digidx == -1 we are out of digits */
151 if (digidx == -1) {
152 break;
153 }
154 /* read next digit and reset the bitcnt */
155 buf = X->dp[digidx--];
156 bitcnt = (int) DIGIT_BIT;
157 }
158
159 /* grab the next msb from the exponent */
160 y = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1;
161 buf <<= (mp_digit)1;
162
163 /* if the bit is zero and mode == 0 then we ignore it
164 * These represent the leading zero bits before the first 1 bit
165 * in the exponent. Technically this opt is not required but it
166 * does lower the # of trivial squaring/reductions used
167 */
168 if (mode == 0 && y == 0) {
169 continue;
170 }
171
172 /* if the bit is zero and mode == 1 then we square */
173 if (mode == 1 && y == 0) {
174 if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
175 goto LBL_RES;
176 }
177 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
178 goto LBL_RES;
179 }
180 continue;
181 }
182
183 /* else we add it to the window */
184 bitbuf |= (y << (winsize - ++bitcpy));
185 mode = 2;
186
187 if (bitcpy == winsize) {
188 /* ok window is filled so square as required and multiply */
189 /* square first */
190 for (x = 0; x < winsize; x++) {
191 if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
192 goto LBL_RES;
193 }
194 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
195 goto LBL_RES;
196 }
197 }
198
199 /* then multiply */
200 if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
201 goto LBL_RES;
202 }
203 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
204 goto LBL_RES;
205 }
206
207 /* empty window and reset */
208 bitcpy = 0;
209 bitbuf = 0;
210 mode = 1;
211 }
212 }
213
214 /* if bits remain then square/multiply */
215 if (mode == 2 && bitcpy > 0) {
216 /* square then multiply if the bit is set */
217 for (x = 0; x < bitcpy; x++) {
218 if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
219 goto LBL_RES;
220 }
221 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
222 goto LBL_RES;
223 }
224
225 bitbuf <<= 1;
226 if ((bitbuf & (1 << winsize)) != 0) {
227 /* then multiply */
228 if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
229 goto LBL_RES;
230 }
231 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
232 goto LBL_RES;
233 }
234 }
235 }
236 }
237
238 mp_exch (&res, Y);
239 err = MP_OKAY;
240 LBL_RES:mp_clear (&res);
241 LBL_MU:mp_clear (&mu);
242 LBL_M:
243 mp_clear(&M[1]);
244 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
245 mp_clear (&M[x]);
246 }
247 return err;
248 }
249 #endif